Nothing Special   »   [go: up one dir, main page]
Next Article in Journal
Acknowledgement to Reviewers of Computation in 2015
Previous Article in Journal
Modeling Groundwater Flow in Heterogeneous Porous Media with YAGMod
You seem to have javascript disabled. Please note that many of the page functionalities won't work as expected without javascript enabled.
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
Computation
Volume 4
Issue 1
10.3390/computation4010003
computation-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

A Test of Various Partial Atomic Charge Models for Computations on Diheteroaryl Ketones and Thioketones

by
Piotr Matczak
Department of Theoretical and Structural Chemistry, Faculty of Chemistry, University of Łódź, Pomorska 163/165, Lodz 90-236, Poland
Computation 2016, 4(1), 3; https://doi.org/10.3390/computation4010003
Submission received: 11 December 2015 / Revised: 2 January 2016 / Accepted: 12 January 2016 / Published: 19 January 2016
(This article belongs to the Section Computational Chemistry)
Browse Figures
Figure 1
<p>The structural skeletal formulas of <b>1a</b>–<b>5a</b> and <b>1b</b>–<b>5b</b>.</p> "> Figure 2
<p>Three conformations of the investigated diheteroaryl ketones and thioketones. The B3LYP/def2-QZVPP-optimized conformers of <b>2a</b> are shown as an example.</p> "> Figure 3
<p>MSE and RMSE (in Debyes) in the μ values approximated by 10 models of partial atomic charge for 30 conformers of <b>1a</b>–<b>5a</b> and <b>1b</b>–<b>5b</b>. The errors are calculated relative to the reference results obtained from the full density at the HF/def2-QZVPP, B3LYP/def2-QZVPP and MP2/def2-QZVPP levels of theory.</p> "> Figure 4
<p>MSE and RMSE (in Debyes) in the μ values approximated by 10 models of partial atomic charge for 30 conformers of <b>1a</b>–<b>5a</b> and <b>1b</b>–<b>5b</b>. The errors are calculated relative to the reference results obtained from the full density at the B3LYP/def2-QZVPP level of theory.</p> ">
Versions Notes

Abstract

:
The effective use of partial atomic charge models is essential for such purposes in molecular computations as a simplified representation of global charge distribution in a molecule and predicting its conformational behavior. In this work, ten of the most popular models of partial atomic charge are taken into consideration, and these models operate on the molecular wave functions/electron densities of five diheteroaryl ketones and their thiocarbonyl analogs. The ten models are tested in order to assess their usefulness in achieving the aforementioned purposes for the compounds in title. Therefore, the following criteria are used in the test: (1) how accurately these models reproduce the molecular dipole moments of the conformers of the investigated compounds; (2) whether these models are able to correctly determine the preferred conformer as well as the ordering of higher-energy conformers for each compound. The results of the test indicate that the Merz-Kollman-Singh (MKS) and Hu-Lu-Yang (HLY) models approximate the magnitude of the molecular dipole moments with the greatest accuracy. The natural partial atomic charges perform best in determining the conformational behavior of the investigated compounds. These findings may constitute important support for the effective computations of electrostatic effects occurring within and between the molecules of the compounds in question as well as similar compounds.

1. Introduction

Partial atomic charges are useful descriptors for interpreting the results of quantum chemical calculations in a chemically intuitive fashion. It comes down to the description of the electron charge distribution within molecules through assigning a partial charge to each atom of the molecules. Because partial atomic charges are not physical observables and there is no strict quantum mechanical definition of them, many different models have been proposed to extract partial atomic charges from the molecular charge distribution. In principle, the models of partial atomic charges may be classified into three groups [1]. The first one covers models based on partitioning the molecular wave function into atomic contributions in terms of basis functions used to construct this wave function. The Mulliken population analysis [2] and the natural population analysis (NPA) [3,4] are typical examples of such models. The Mulliken population analysis is probably the best known of all models of partial atomic charge. Due to its simplicity, this model is computationally very attractive. However, its results tend to vary with the basis set employed and they are unrealistic in some cases [5]. These drawbacks partially arise from the fact that the Mulliken population analysis utilizes a nonorthogonal basis set. This problem is overcome by the NPA in which orthonormal natural atomic functions are constructed. The distributed multipole analysis of Gaussian wave functions (GDMA) [6] is also ranked among the models of the first group. In this analysis, the distribution of the molecular electron density is represented by multipoles located on the individual atoms of a molecule. The partitioning of the molecular electron density between the atoms is carried out in the basis-function space because the multipoles are obtained by using the products of Gaussian basis functions. The second group of partial atomic charge models comprises models that partition the molecular electron density into atomic domains in the physical space. These atomic domains are defined by means of Bader’s atoms-in-molecules (AIM) topological analysis of the molecular electron density [7] or using the so-called “promolecular” density proposed originally by Hirshfeld [8]. To be more specific, the AIM division of a molecule into atoms is based on the zero-flux surfaces of the molecular electron density, and then this density is integrated over the resulting atomic domains to obtain their partial charges. In the Hirshfeld model the partial charge of each atom is calculated by assuming that the electron density at each point is shared between the surrounding atoms in direct proportion to their free-atom electron densities at the corresponding distances from the nuclei. The original partial atomic charges obtained from the Hirshfeld population analysis can be improved through parametrization to reproduce accurately a particular molecular property. For instance, the recently developed charge model 5 (CM5) utilizes a set of parameters derived by fitting to reference values of gas-phase dipole moments [9]. As for the models belonging to the third group, their partial atomic charges are based on the reproduction of the molecular electrostatic potential. Various schemes that fit partial atomic charges to the molecular electrostatic potential are reported, e.g., the Merz-Kollman-Singh (MKS) scheme [10,11], charges from electrostatic potentials (CHELP) [12], charges from electrostatic potentials using a grid (CHELPG) [13] and the Hu-Lu-Yang (HLY) scheme [14]. The MKS, CHELP and CHELPG schemes differ in the selection of points surrounding a molecule. Such points are used for the calculation of the molecular electrostatic potential in the fitting procedure. The MKS scheme employs points located at four shells at a distance of 1.4, 1.6, 1.8 and 2.0 times the van der Waals radii of the atoms constituting a molecule. The CHELP scheme samples the molecular electrostatic potential at points at a distance of 2.5, 3.5, 4.5, 5.5 and 6.5 Å from the van der Waals surface. In the CHELPG scheme the points are selected from a regularly spaced grid, between 0 and 2.8 Å from the van der Waals surface. The HLY scheme introduces the so-called object function in the entire molecular volume space instead of the points around the molecule. The application of the object function improves the numerical stability of fitting results.
Partial atomic charges find application in the molecular modeling of several chemically relevant areas, such as explaining structural and reactivity differences between various molecules and their conformers [15,16,17], investigating charge transfers within a single molecule and between several molecules [18,19,20], and predicting pKa variations in series of molecules [21,22]. Since the definition of partial atomic charges is not strict, various models of partial atomic charges can, however, differ significantly in the reliability of their predictions. Therefore, an evaluation of the performance of a partial atomic charge model for a given problem is necessary before using it for making reliable predictions [23,24,25,26,27,28].
In this work, various models of partial atomic charge are tested to establish their efficiency for computations on diheteroaryl ketones and thioketones. A series of five diheteroaryl ketones and their thiocarbonyl analogs is considered. These ketones and thioketones are formed by the disubstitution of formaldehyde and thioformaldehyde with a five-membered heterocyclic group. The series includes di(furan-2-yl)methanone (1a), di(thiophen-2-yl)methanone (2a), di(selenophen-2-yl)methanone (3a), di(1H-pyrrol-2-yl)methanone (4a), di(1-methylpyrrol-2-yl)methanone (5a) and their thiocarbonyl analogs (1b–5b). The structural skeletal formulas of 1a–5a and 1b–5b are shown in Figure 1. The thioketones 1b–5b have recently been synthesized by the O/S exchange in the corresponding ketones by treatment with Lawesson’s reagent [29]. In the first stage of the present work, various models of partial atomic charge are tested in terms of their ability to reproduce the molecular dipole moments of 1a–5a and 1b–5b. Next, the efficiency of the models for predicting the conformational behavior of 1a–5a and 1b–5b is studied. The prediction of the conformational behavior will be restricted here to the determination of a preferred conformer and the sequence of higher-energy conformers.
Computation 04 00003 g001 1024
Figure 1. The structural skeletal formulas of 1a–5a and 1b–5b.
Figure 1. The structural skeletal formulas of 1a–5a and 1b–5b.
Computation 04 00003 g001

2. Computational Details

The molecular structures of 1a–5a and 1b–5b are taken from our previous works [30,31,32] in which Becke’s three-parameter hybrid exchange functional combined with the correlation functional of Lee, Yang and Parr (B3LYP) [33,34,35] and the def2-QZVPP basis set [36] were used to optimize the geometries of the isolated molecules of 1a–5a and 1b–5b. It was also established there that for each of the compounds its three conformations could be formed by the rotation of heteroaryl substituents about the single C-C bonds linking these substituents with the C atom of carbonyl/thiocarbonyl group. The resulting conformers of 1a–5a and 1b–5b are denoted by the prefixes cc, ct and tt, indicating the spatial arrangement of ring heteroatoms with respect to the O/S atom of carbonyl/thiocarbonyl group (see Figure 2).
Computation 04 00003 g002 1024
Figure 2. Three conformations of the investigated diheteroaryl ketones and thioketones. The B3LYP/def2-QZVPP-optimized conformers of 2a are shown as an example.
Figure 2. Three conformations of the investigated diheteroaryl ketones and thioketones. The B3LYP/def2-QZVPP-optimized conformers of 2a are shown as an example.
Computation 04 00003 g002
For the cc-, ct- and tt-conformers of 1a–5a and 1b–5b, their molecular wave functions/electron densities are calculated at three levels of theory, that is, HF/def2-QZVPP [36,37,38], B3LYP/def2-QZVPP [33,34,35,36] and MP2/def2-QZVPP [36,39]. These molecular wave functions/electron densities are the starting point for deriving partial atomic charges from various models belonging to the three groups mentioned in the introduction. The partial atomic charges on all atoms of each conformer are determined by means of ten models, namely Mulliken, NPA, GDMA, AIM, Hirshfeld, CM5, MKS, CHELP, CHELPG and HLY. Then, the magnitude of the dipole moment μ of the conformer is approximated by the following formula:
μ = ∑ Q = x , y , z ( ∑ k = 1 N q k Q k ) 2
where the index k runs over all atoms of the conformer having N atoms, Q denotes the Cartesian coordinates of the corresponding atom and q is the partial charge of the atom. In order to assess the accuracy of the dipole moments calculated using Equation (1), their values need to be compared with the corresponding reference results, preferably obtained from experimental measurements. However, such reference results are not available for individual gas-phase conformers of 1a–5a and 1b–5b. Therefore, the dipole moments calculated by the ten models of partial atomic charge are compared with the corresponding dipole moments obtained from the regular quantum chemical calculations involving the full molecular electron density (to be strict, these dipole moments are calculated as an expectation value of the appropriate quantum operator, which is well defined for the dipole moment).
The calculated partial atomic charges are also used for a rough estimation of the magnitude of electrostatic effects occurring in the conformers of 1a–5a and 1b–5b. On this basis, the conformational behavior of these compounds can be predicted. Interaction energies in the pairs of partial atomic charges (EC(i,j)) are calculated using classical Coulomb’s law. These interaction energies are summed up over all pairs of partial atomic charges within each conformer (Σi>j EC(i,j)) to roughly estimate the energy Eelst associated with electrostatic effects occurring in the conformer. The conformers of each compound can be ordered with respect to their Eelst values. In the resulting sequence, the lower (that is, the more negative) the value of Eelst is obtained for a conformer, the more stable the conformer is. Such a procedure taking only Eelst into consideration assumes that the electrostatic effects play an important role in governing the conformational behavior of investigated compounds. These effects, indeed, contribute significantly to the ordering of conformers for diheteroaryl ketones and thioketones [30,32].
The GDMA 2.2.09 [6] and AIMAll 14.06.21 [40] programs have been used to calculate partial atomic charges derived from the GDMA and AIM models, respectively. All the remaining calculations have been carried out with the Gaussian 09 D.01 program [41]. In this program, the fitting of MKS, CHELP, CHELPG and HLY partial atomic charges to reproduce the molecular electrostatic potential has involved no additional constraint of reproducing the molecular dipole moment.

3. Results and Discussion

Since the molecular dipole moment is a primary quantity providing an essential insight into the distribution of electron charge within a molecule, a reasonable test for any model of partial atomic charge is to inspect the performance of the model in predicting such a quantity. Therefore, we start testing ten models of partial atomic charge by establishing their ability to reproduce the magnitude of the dipole moments for the conformers of 1a–5a and 1b–5b. Figure 3 shows the mean signed error (MSE) and root mean square error (RMSE) in the μ values approximated by the partial atomic charges derived from each model with respect to the reference results obtained from the full-density calculations. Additionally, the values of MSE and RMSE are determined for three levels of theory. When comparing here the approximated values of μ with the corresponding reference values, both the former and the latter are obtained from the molecular wave functions/electron densities generated at the same level of theory. The approximated and reference μ values used for the calculations of the MSE and RMSE presented in Figure 3 can be found in Tables S1–S6 in Supplementary Materials. The MSE provides information about systematic errors occurring in the approximated values of μ, whereas the RMSE allows us to rank the accuracy of individual models for reproducing the reference values of μ.
It is evident from what the lower plot in Figure 3 shows that the μ values obtained from the MKS and HLY partial atomic charges reproduce best the μ values calculated from the full density. Among the models that are not based on electrostatic potential fitting, the CM5 model offers the highest accuracy of the approximated μ values. The μ values approximated by the AIM model lead to the largest RMSE values. The accuracy of four models that are based on the molecular electrostatic potential is not very differentiated: the RMSE values of these models fall in the range from 0.03 D (for the HLY model) to 0.48 D (for the CHELP model). As indicated by the values of the MSE in the upper plot in Figure 3, the AIM and GDMA models overestimate the magnitude of μ significantly, whereas the Mulliken, Hirshfeld, CM5 and CHELP models tend to underestimate the values of μ. The differences in the MSE and RMSE values obtained from various levels of theory are usually relatively small, and thus all the above-mentioned findings are common to the HF, B3LYP and MP2 levels.
Computation 04 00003 g003 1024
Figure 3. MSE and RMSE (in Debyes) in the μ values approximated by 10 models of partial atomic charge for 30 conformers of 1a–5a and 1b–5b. The errors are calculated relative to the reference results obtained from the full density at the HF/def2-QZVPP, B3LYP/def2-QZVPP and MP2/def2-QZVPP levels of theory.
Figure 3. MSE and RMSE (in Debyes) in the μ values approximated by 10 models of partial atomic charge for 30 conformers of 1a–5a and 1b–5b. The errors are calculated relative to the reference results obtained from the full density at the HF/def2-QZVPP, B3LYP/def2-QZVPP and MP2/def2-QZVPP levels of theory.
Computation 04 00003 g003
Our findings should now be collated with previous observations on the performance of various models of partial atomic charge. It is known that the NPA partial atomic charges overestimate the dipole moments of polypeptides [42] and in the present work this observation is validated for smaller molecular systems (see the positive values of the MSE for the NPA model in Figure 3). The overestimation of the μ values predicted by the NPA model seems to be associated with the fact that this model generally overestimates the amplitude of the electrostatic potential, as it was reported for a large set of 500 simple organic molecules [43]. In another work taking a large set of simple molecules into account [9], it was found that the molecular dipole moment is usually underestimated when it is approximated by the Hirshfeld model. It is valid in our case, as evidenced by the negative MSE values of the Hirshfeld model in Figure 3. Both the MSE and RMSE values in Figure 3 indicate that the CM5 model outperforms the Hirshfeld one, which is also in accord with previous reports [9,44]. Furthermore, the CM5 model turns out to be slightly more accurate than the Mulliken one, and this observation is shared by the results of dipole moment calculations using the MG3S basis set [9]. It should, however, be stressed here that the CM5 model is essentially independent of a basis set [9], while the basis set dependence is the well-known weakness of the Mulliken model [1,4]. The poor performance of the AIM model results from its excessive partial atomic charges. The severe overestimation of dipole moment magnitude is typical of this model and it is observed even for the water molecule [45]. As noted in the previous paragraph, the MKS and HLY models provide excellent results for the μ values of the conformers of 1a–5a and 1b–5b. The CHELP and CHELPG models demonstrate larger RMSE values but both are still fairly successful in reproducing the μ values from the full density. This seems to be in line with what was previously reported for the MKS and CHELPG models [46]. The former turned out to be superior to the latter for the representation of dipole moments in ionic liquids. One of the reasons for the good performance of the MKS, CHELP, CHELPG and HLY models in reproducing μ for the conformers of 1a–5a and 1b–5b is that the molecular shape of these conformers is not very complex (the heteroaryl substituents are planar although they are most often not coplanar with one another) and almost all their atoms are near their molecular van der Waals surfaces. This practically eliminates the occurrence of the so-called “buried atoms” for which the electrostatic potential fitting is inaccurate, and thus, the resulting partial atomic charges are poorly determined [47].
Aside from the statistical comparison of the approximated μ values with the reference results from the full densities computed at the HF/def2-QZVPP, B3LYP/def2-QZVPP or MP2/def2-QZVPP level of theory, another comparison of the calculated μ values is meaningful. Such a comparison utilizes only a single set of reference μ values, irrespective of the level of theory which produces the molecular wave functions/electron densities for the models of partial atomic charge. This comparison allows us to directly examine whether the level of theory affects the performance of the ten models in predicting the values of μ. Because the range of the quantum chemical methods used in this work includes both the simplest ab initio method (that is, HF) and two more advanced methods (that is, B3LYP and MP2), the effect of electron correlation on the performance of the models can be studied. The B3LYP/def2-QZVPP level of theory is selected to provide the reference full-density values of μ for the conformers of 1a–5a and 1b–5b. Of three levels considered in this work, the μ values obtained from the B3LYP/def2-QZVPP full density turned out to be closest to experiment for a test set of simple molecules being the building blocks of 1a–5a and 1b–5b (see section S2 in Supplementary Materials). On this basis, the μ values calculated using the B3LYP/def2-QZVPP full density [32] are also assumed to be the most realistic for the conformers of 1a–5a and 1b–5b. These values are now used as the only reference results for the calculations of the MSE and RMSE in the μ values approximated by the ten models that, in turn, operate on the HF/def2-QZVPP, B3LYP/def2-QZVPP and MP2/def2-QZVPP wave functions/electron densities. The calculated MSE and RMSE values are presented graphically in Figure 4. For the B3LYP method, its results shown in this figure are obviously identical to those depicted in Figure 3.
Results shown in Figure 4 support our previous findings that such electrostatic potential-based models as MKS, CHELP, CHELPG and HLY are able to predict the values of μ with great accuracy. Out of these four models, it is hard to select a single model that performs best at all three levels of theory. For the molecular wave functions calculated by HF/def2-QZVPP, the CHELP model leads to the lowest RMSE value. For the molecular wave functions computed at the B3LYP/def2-QZVPP and MP2/def2-QZVPP levels, the most accurate values of μ are predicted by the HLY and MKS models, respectively. The positive values of the MSE yielded by the MKS, CHELP, CHELPG and HLY models operating on the HF/def2-QZVPP wave functions illustrate the systematic overestimation of the approximated μ values. This is in line with the previous finding made for the MKS model [47], and now it is additionally extended to the CHELP, CHELPG and HLY models. All four models consistently change their behavior while employing the MP2/def2-QZVPP wave functions (MSE < 0). Differences in the three RMSE values of each model are usually relatively small and in such cases they do not exceed 0.5 D. This indicates that, in general, the effect of the quantum chemical method applied for obtaining the molecular wave functions/electron densities is rather minor. However, at least one exception to this general conclusion can be seen in Figure 4. The AIM model combined with the HF/def2-QZVPP electron densities performs particularly badly when compared with its RMSE values obtained from both B3LYP/def2-QZVPP and MP2/def2-QZVPP electron densities. From this perspective, the AIM model seems to be particularly sensitive to the inclusion of electron correlation effects into the molecular electron density. On the other hand, it should be stressed that the performance of the AIM model still remains poor for the B3LYP/def2-QZVPP and MP2/def2-QZVPP electron densities.
Computation 04 00003 g004 1024
Figure 4. MSE and RMSE (in Debyes) in the μ values approximated by 10 models of partial atomic charge for 30 conformers of 1a–5a and 1b–5b. The errors are calculated relative to the reference results obtained from the full density at the B3LYP/def2-QZVPP level of theory.
Figure 4. MSE and RMSE (in Debyes) in the μ values approximated by 10 models of partial atomic charge for 30 conformers of 1a–5a and 1b–5b. The errors are calculated relative to the reference results obtained from the full density at the B3LYP/def2-QZVPP level of theory.
Computation 04 00003 g004
A brief search of the literature reveals several previous mentions of the impact of electron correlation on partial atomic charges and resulting μ values. It was reported that this impact is small [48] or even minimal [46]. Our results are in agreement with these findings. The usual effect of electron correlation on the distribution of molecular electron density is to deplete electron density from the centers of bonds and increase it in shells around the atomic nuclei and at periphery of molecules [49]. However, the influence of such a depletion on molecular dipole moments is often smaller than might be expected because the charge reorganization happens around each atom and can result in only a small net change in polarization of whole molecules [48]. In the case of the AIM model, a significant insensitivity of its many parameters (also those obtained through the integration of electron density) to the addition of electron correlation and the change in basis set was observed [50,51]. On the other hand, it was detected for the water molecule that, within the framework of the AIM model, the electron density produced by the HF method leads to a slightly enhanced dipole moment compared to that obtained from the MP2 electron density [45]. This observation is in agreement with our results that, indeed, show an increase in the RMSE values for the AIM model while going from the correlation-corrected electron densities to the HF/def2-QZVPP densities.
The next stage of the present work is to establish the usability of various partial atomic charge models to predict qualitatively the conformational behavior of 1a–5a and 1b–5b. The Eelst energy calculated using the partial atomic charges determined for each conformer is considered here to be a simple measure of electrostatic effects stabilizing the conformer. It is convenient to express the Eelst values obtained for the conformers of each compound with respect to the Eelst energy of the conformer that is most favorable (in other words, with respect to the conformer possessing the lowest Eelst energy). The resulting relative electrostatic energy ΔEelst is equal to zero for the preferred conformer while it is positive for less stable conformers. A part of the ΔEelst values characterizing the cc-, ct-, and tt-conformers of 1a–5a and 1b–5b is given in Table 1 and Table 2. The tabulated values have been calculated using selected models that have operated on the molecular wave functions/electron densities computed at the B3LYP/def2-QZVPP level of theory (a complete set of results obtained from all ten models of partial atomic charge, as well as at the HF/def2-QZVPP, B3LYP/def2-QZVPP and MP2/def2-QZVPP levels of theory can be found in Tables S8–S13 in Supplementary Materials). Additionally, the relative electron energies ΔE and full-density dipole moments μ obtained from the previous regular calculations at the B3LYP/def2-QZVPP level [30,32] are also presented in Table 1 and Table 2. The values of ΔE allow us to establish the reference orderings of conformers for 1a–5a and 1b–5b.
Table
Table 1. Relative electron energies (ΔE in kcal/mol), dipole moments obtained from the full density (μ in Debyes) and relative electrostatic energies (ΔEelst in kcal/mol) for 1a–5a in their three conformations. All the results are calculated at the B3LYP/def2-QZVPP level of theory.
Table 1. Relative electron energies (ΔE in kcal/mol), dipole moments obtained from the full density (μ in Debyes) and relative electrostatic energies (ΔEelst in kcal/mol) for 1a–5a in their three conformations. All the results are calculated at the B3LYP/def2-QZVPP level of theory.
ConformerΔE aμ bΔEelst
MullikenNPAAIMHirshfeldMKSHLY
cc-1a2.144.7313.788.8035.421.810.070.00
ct-1a0.003.914.130.000.951.100.4719.65
tt-1a0.342.910.000.670.000.000.0025.15
cc-2a0.004.030.930.000.000.000.000.00
ct-2a0.773.370.305.715.290.6844.4764.19
tt-2a1.852.690.0010.9910.581.4750.2869.59
cc-3a0.003.650.000.000.000.000.000.00
ct-3a1.413.211.017.6010.110.1553.4668.89
tt-3a2.932.811.4216.4422.960.0569.2587.14
cc-4a0.000.410.000.005.380.000.000.00
ct-4a3.823.3110.997.790.001.756.791.51
tt-4a8.995.319.7711.968.283.3426.9225.70
cc-5a0.000.130.000.003.350.000.000.00
ct-5a5.373.558.001.813.300.5036.1944.15
tt-5a9.635.2928.215.540.001.3699.21106.26
a Values taken from [30]; b Values taken from [32].
The tabulated values of ΔEelst show that the energy of electrostatic interactions between partial atomic charges is able to give us some indication of the preferred conformations for the investigated compounds, especially for those whose preferred conformation is characterized by either the largest or the smallest values of μ. The most energetically stable conformers of 1a and 1b possess molecular dipole moments whose values lie in the middle between the values obtained for the other two conformers. In such case the Eelst energy usually turns out to be insufficient to identify the preferred conformation. Although the preference of the ct-conformation can be inferred from ΔEelst calculated using the NPA and Mulliken partial atomic charges, none of the two sets of partial atomic charges leads to ct-conformation preference simultaneously for 1a and 1b. Nevertheless, the NPA model is recognized to be most successful in predicting the preferred conformers of 1a–5a and 1b–5b in terms of Eelst. The Mulliken, AIM and Hirshfeld models lead to a slightly greater number of incorrect indications of preferred conformers than the NPA model does. Furthermore, this model always performs best, irrespective of the level of theory used to generate the molecular wave functions. The electrostatic potential-based models are generally less accurate in identifying the preferred conformers of 1a–5a and 1b–5b than the models belonging to the two remaining classes. Of the electrostatic potential-based models, only MKS and HLY are able to correctly indicate the preferred conformers for more than half of the investigated compounds.
Table
Table 2. Relative electron energies (ΔE in kcal/mol), dipole moments obtained from the full density (μ in Debyes) and relative electrostatic energies (ΔEelst in kcal/mol) for 1b–5b in their three conformations. All the results are calculated at the B3LYP/def2-QZVPP level of theory.
Table 2. Relative electron energies (ΔE in kcal/mol), dipole moments obtained from the full density (μ in Debyes) and relative electrostatic energies (ΔEelst in kcal/mol) for 1b–5b in their three conformations. All the results are calculated at the B3LYP/def2-QZVPP level of theory.
ConformerΔE aμ bΔEelst
MullikenNPAAIMHirshfeldMKSHLY
cc-1b1.464.697.790.000.000.8051.0242.63
ct-1b0.004.090.001.352.220.6210.3511.94
tt-1b0.313.350.7510.9036.070.000.000.00
cc-2b0.004.130.000.000.000.000.000.00
ct-2b0.993.621.022.842.090.4837.2549.47
tt-2b2.223.081.236.274.481.1535.8251.96
cc-3b0.003.800.000.000.000.250.000.00
ct-3b1.443.462.093.291.510.2756.3257.83
tt-3b3.003.173.058.545.660.0065.3167.72
cc-4b0.001.540.640.000.000.0027.9431.98
ct-4b3.923.965.861.5911.121.2237.8241.82
tt-4b9.415.660.000.0138.501.950.000.00
cc-5b0.001.230.000.000.000.00180.49180.68
ct-5b3.054.4410.212.660.010.91117.01113.67
tt-5b5.045.6217.642.617.421.600.000.00
a Values taken from [30]; b Values taken from [32].
Aside from the application of Eelst to the identification of the most stable conformation, this quantity allows us to order conformers relative to their Eelst values and the resulting orderings of conformers for 1a–5a and 1b–5b can be compared with the orderings based on ΔE. The evident relationship between μ and ΔE for 2a–5a and their thiocarbonyl counterparts suggests that electrostatic effects are responsible to a certain extent for the ordering of individual conformers relative to their energetic stability. We focus here only on the successful reproduction of conformer orderings and not on the quantitative correlation between individual non-zero ΔEelst and ΔE values. The values of ΔEelst should not be compared directly with ΔE because the former are merely a crude approximation of intramolecular electrostatic effects. The electrostatic effects are undoubtedly important but not the sole factor affecting the stability of the investigated conformers. Therefore, the values of ΔEelst are usually far from the corresponding ΔE values. The values of ΔEelst obviously inherit the deficiencies of the applied model of partial atomic charge. The ΔEelst values calculated using the AIM partial atomic charges are usually large for 1a–3a and 1b–3b, whereas the Hirshfeld partial atomic charges lead to very small ΔEelst values for the majority of the investigated compounds. This is due to the fact that AIM partial atomic charges generally tend to adopt large absolute values [23], while the Hirshfeld model assigns partial atomic charges that are very small in magnitude [52].
Turning our attention to the predicted conformer orderings, we can see that none of the considered models of partial atomic charge provides the correct orderings of conformers for all ten compounds. Among the ten models, the NPA one appears to be most suitable for illustrating tendencies in electrostatic effects (in terms of Eelst) in the investigated compounds, although some inconsistencies with the orderings predicted by ΔE are found. The inconsistencies in the conformer orderings yielded by the NPA model occur mostly for diheteroaryl thioketones. The AIM and Hirshfeld models give less satisfactory results than those of the NPA model. On the other hand, they reproduce the reference orderings of conformers for a greater number of the investigated compounds than the HLY model does. Of the electrostatic potential-based models, HLY produces the conformer orderings that fit best to the corresponding reference orderings indicated by ΔE. All the aforementioned findings are common to the three levels of theory used to generate the wave functions/electron densities of the conformers. This is additionally supported by the results presented in Table 3. This table shows the percentage similarity of the conformer sequences deduced from ΔEelst to the reference conformer orderings determined in terms of ΔE. The percentage similarity has been obtained through comparing the conformer orderings of all investigated compounds.
Table
Table 3. Percentage similarity of the conformer sequences predicted by ΔEelst to the reference sequences determined using ΔE. The values of ΔEelst are calculated by 10 partial atomic charge models operating on the wave functions/electron densities calculated at three levels of theory.
Table 3. Percentage similarity of the conformer sequences predicted by ΔEelst to the reference sequences determined using ΔE. The values of ΔEelst are calculated by 10 partial atomic charge models operating on the wave functions/electron densities calculated at three levels of theory.
ModelLevel of Theory
HFB3LYPMP2
Mulliken507050
NPA837383
GDMA575357
AIM677080
Hirshfeld737070
CM5374347
MKS606050
CHELP233320
CHELPG574353
HLY606763

4. Conclusions

In this work, ten of the most popular models of partial atomic charge have been considered, and these models have operated on the molecular wave functions/electron densities of the conformers of 1a–5a and 1b–5b. These wave functions/electron densities have been calculated at three different levels of theory. The ten models were tested in order to assess their usefulness in performing effective computations on diheteroaryl ketones and thioketones. More specifically, our test assesses the models’ abilities (1) to approximate the magnitude of μ for the conformers of 1a–5a and 1b–5b, and (2) to correctly determine the conformers’ orderings through the estimation of the electrostatic interaction between partial atomic charges within the conformers.
The results of the test indicate that the simplified representation of the magnitude of μ by partial atomic charges is most effective when the MKS and HLY models are used to produce the partial atomic charges within the conformers of 1a–5a and 1b–5b. These models are able to reproduce very accurately the reference μ values obtained from the full density, and they perform well for the molecular wave functions calculated at all three levels of theory. Among the models that are not based directly on the molecular electrostatic potential, the CM5 model offers a reasonable accuracy in approximating the values of μ. From the subsequent results of our test it can be concluded that the most successful estimation of the intramolecular electrostatic effects governing the conformational behavior of 1a–5a and 1b–5b is provided by the partial atomic charges derived from the NPA model. Besides the designation of the most successful model for the determination of the conformational behavior, this part of the test also shows that the simple approach utilizing the calculation of Eelst by means of NPA partial atomic charges is a surprisingly effective yet still qualitative tool to anticipate the energetic orderings of the conformers of 1a–5a and 1b–5b. It also implies an important role of intramolecular electrostatic effects in determining the conformational behavior of the investigated compounds.
The above-mentioned results of our test may be essential for establishing the effective approximations of molecular dipole moments and intramolecular electrostatic effects for future computations on the molecules of not only the compounds in question but also similar compounds. Furthermore, these results may have implications for the tuning of force fields used in the classical molecular dynamics simulations of diheteroaryl ketones and thioketones. The accurate reproduction of molecular dipole moments by partial atomic charges is an important factor contributing to the reliable determination of intermolecular interactions in such simulations.

Supplementary Files

Supplementary File 1

Acknowledgments

This work was partially supported by PL-Grid Infrastructure.

Conflicts of Interest

The author declares no conflict of interest.

References

  1. Jensen, F. Introduction to Computational Chemistry, 2nd ed.; John Wiley & Sons, Inc.: Chichester, UK, 2007. [Google Scholar]
  2. Mulliken, R.S. Electronic population analysis on LCAO-MO molecular wave functions. J. Chem. Phys. 1955, 23, 1833–1840. [Google Scholar] [CrossRef]
  3. Reed, A.E.; Curtiss, L.A.; Weinhold, F. Intermolecular interactions from a natural bond orbital, donor-acceptor viewpoint. Chem. Rev. 1988, 88, 899–926. [Google Scholar] [CrossRef]
  4. Reed, A.E.; Weinstock, R.B.; Weinhold, F. Natural population analysis. J. Chem. Phys. 1985, 83, 735–746. [Google Scholar] [CrossRef]
  5. Thompson, J.D.; Xidos, J.D.; Sonbuchner, T.M.; Cramer, C.J.; Truhlar, D.G. More reliable partial atomic charges when using diffuse basis sets. PhysChemComm 2002, 5, 117–134. [Google Scholar] [CrossRef]
  6. Stone, A.J. Distributed multipole analysis: Stability for large basis sets. J. Chem. Theory Comput. 2005, 1, 1128–1132. [Google Scholar] [CrossRef] [PubMed]
  7. Bader, R.F.W. Atoms in Molecules: A Quantum Theory; Clarendon: Oxford, UK, 1990. [Google Scholar]
  8. Hirshfeld, F.L. Bonded-atom fragments for describing molecular charge densities. Theor. Chem. Acta 1977, 44, 129–138. [Google Scholar] [CrossRef]
  9. Marenich, A.V.; Jerome, S.V.; Cramer, C.J.; Truhlar, D.G. Charge model 5: An extension of Hirshfeld population analysis for the accurate description of molecular interactions in gaseous and condensed phases. J. Chem. Theory Comput. 2012, 8, 527–541. [Google Scholar] [CrossRef] [PubMed]
  10. Singh, U.C.; Kollman, P.A. An approach to computing electrostatic charges for molecules. J. Comput. Chem. 1984, 5, 129–145. [Google Scholar] [CrossRef]
  11. Besler, B.H.; Merz, K.M., Jr.; Kollman, P.A. Atomic charges derived from semiempirical methods. J. Comput. Chem. 1990, 11, 431–439. [Google Scholar] [CrossRef]
  12. Chirlian, L.E.; Francl, M.M. Atomic charges derived from electrostatic potentials—A detailed study. J. Comput. Chem. 1987, 8, 894–905. [Google Scholar] [CrossRef]
  13. Breneman, C.M.; Wiberg, K.B. Determining atom-centered monopoles from molecular electrostatic potentials—The need for high sampling density in formamide conformational analysis. J. Comput. Chem. 1990, 11, 361–373. [Google Scholar] [CrossRef]
  14. Hu, H.; Lu, Z.; Yang, W. Fitting molecular electrostatic potentials from quantum mechanical calculations. J. Chem. Theory Comput. 2007, 3, 1004–1013. [Google Scholar] [CrossRef] [PubMed]
  15. Rostkowski, M.; Paneth, P. Charge localization in monothiophosphate monoanions. Pol. J. Chem. 2007, 81, 711–720. [Google Scholar]
  16. Balasubramanian, G.; Schulte, J.; Müller-Plathe, F.; Böhm, M.C. Structural and thermochemical properties of a photoresponsive spiropyran and merocyanine pair: Basis set and solvent dependence in density functional predictions. Chem. Phys. Lett. 2012, 554, 60–66. [Google Scholar] [CrossRef]
  17. Cvijetic, I.N.; Vitorovic-Todorovic, M.D.; Juranic, I.O.; Nakarada, D.J.; Milosavljevic, M.D.; Drakulic, B.J. Reactivity of (E)-4-aryl-4-oxo-2-butenoic acid phenylamides with piperidine and benzylamine: Kinetic and theoretical study. Monatsh. Chem. 2014, 145, 1297–1306. [Google Scholar] [CrossRef]
  18. Domagała, M.; Matczak, P.; Palusiak, M. Halogen bond, hydrogen bond and N···C interaction—On interrelation among these three noncovalent interactions. Comput. Theor. Chem. 2012, 998, 26–33. [Google Scholar] [CrossRef]
  19. Domagała, M.; Palusiak, M. The influence of substituent effect on noncovalent interactions in ternary complexes stabilized by hydrogen-bonding and halogen-bonding. Comput. Theor. Chem. 2014, 1027, 173–178. [Google Scholar] [CrossRef]
  20. Guidara, S.; Feki, H.; Abid, Y. Structural, vibrational, NLO, MEP, NBO analysis and DFT calculation of bis 2,5-dimethylanilinium sulfate. J. Mol. Struct. 2015, 1080, 176–187. [Google Scholar] [CrossRef]
  21. Varekova, R.S.; Geidl, S.; Ionescu, C.-M.; Skrehota, O.; Bouchal, T.; Sehnal, D.; Abagyan, R.; Koca, J. Predicting pKa values from EEM atomic charges. J. Cheminform. 2013, 5, 18. [Google Scholar] [CrossRef] [PubMed]
  22. Ugur, I.; Marion, A.; Parant, S.P.; Jensen, J.H.; Monard, G. Rationalization of the pKa values of alcohols and thiols using atomic charge descriptors and its application to the prediction of amino acid pKa’s. J. Chem. Inf. Model. 2014, 54, 2200–2213. [Google Scholar] [CrossRef] [PubMed]
  23. Gross, K.C.; Seybold, P.G.; Hadad, C.M. Comparison of different atomic charge schemes for predicting pKa variations in substituted anilines and phenols. Int. J. Quantum Chem. 2002, 90, 445–458. [Google Scholar] [CrossRef]
  24. Choi, Y.; Kim, H.; Cho, K.W.; Paik, S.R.; Kim, H.-W.; Jeong, K.; Jung, S. Systematic probing of an atomic charge set of sialic acid disaccharides for the rational molecular modeling of avian influenza virus based on molecular dynamics simulations. Carbohydr. Res. 2009, 344, 541–544. [Google Scholar] [CrossRef] [PubMed]
  25. Jacquemin, D.; le Bahers, T.; Adamo, C.; Ciofini, I. What is the ‘‘best’’ atomic charge model to describe through-space charge-transfer excitations? Phys. Chem. Chem. Phys. 2012, 14, 5383–5388. [Google Scholar] [CrossRef] [PubMed]
  26. Jambeck, J.P.M.; Mocci, F.; Lyubartsev, A.P.; Laaksonen, A. Partial atomic charges and their impact on the free energy of solvation. J. Comput. Chem. 2013, 34, 187–197. [Google Scholar] [CrossRef] [PubMed]
  27. Wang, B.; Li, S.L.; Truhlar, D.G. Modeling the partial atomic charges in inorganometallic molecules and solids and charge redistribution in lithium-ion cathodes. J. Chem. Theory Comput. 2014, 10, 5640–5650. [Google Scholar] [CrossRef] [PubMed]
  28. Hamad, S.; Balestra, S.R.G.; Bueno-Perez, R.; Calero, S.; Ruiz-Salvador, A.R. Atomic charges for modeling metal–organic frameworks: Why and how. J. Solid State Chem. 2015, 223, 144–151. [Google Scholar] [CrossRef]
  29. Mlostoń, G.; Urbaniak, K.; Gębicki, K.; Grzelak, P.; Heimgartner, H. Hetaryl thioketones: Synthesis and selected reactions. Heteroat. Chem. 2014, 25, 548–555. [Google Scholar] [CrossRef]
  30. Matczak, P.; Domagała, M.; Domagała, S. Conformers of diheteroaryl ketones and thioketones: A quantum chemical study of their properties and fundamental intramolecular energetic effects. Struct. Chem. 2015. [Google Scholar] [CrossRef]
  31. Matczak, P. Intramolecular C–H···H–C contacts in diheteroaryl ketones and thioketones: A theoretical analysis. Bull. Chem. Soc. Jpn. 2015. [Google Scholar] [CrossRef]
  32. Matczak, P.; Domagała, M. Charge distribution in conformers of diheteroaryl ketones and thioketones. Submitted.
  33. Becke, A.D. Density-functional thermochemistry. III. The role of exact exchange. J. Chem. Phys. 1993, 98, 5648–5642. [Google Scholar] [CrossRef]
  34. Vosko, S.H.; Wilk, L.; Nusair, M. Accurate spin-dependent electron liquid correlation energies for local spin density calculations: A critical analysis. Can. J. Phys. 1980, 58, 1200–1211. [Google Scholar] [CrossRef]
  35. Lee, C.; Yang, W.; Parr, R.G. Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density. Phys. Rev. B 1988, 37, 785–789. [Google Scholar] [CrossRef]
  36. Weigend, F.; Ahlrichs, R. Balanced basis sets of split valence, triple zeta valence and quadruple zeta valence quality for H to Rn: Design and assessment of accuracy. Phys. Chem. Chem. Phys. 2005, 7, 3297–3305. [Google Scholar] [CrossRef] [PubMed]
  37. Hartree, D.R. The wave mechanics of an atom with a non-Coulomb central field. Part I: Theory and methods. Proc. Camb. Philos. Soc. 1928, 24, 89–110. [Google Scholar] [CrossRef]
  38. Fock, V. Näherungsmethode zur lösung des quantenmechanischen mehrkörperproblems. Z. Phys. 1930, 61, 126–148. (In German) [Google Scholar] [CrossRef]
  39. Møller, C.; Plesset, M.S. Note on an approximation treatment for many-electron systems. Phys. Rev. 1934, 46, 618–622. [Google Scholar] [CrossRef]
  40. Keith, T.A. AIMAll (Version 14.06.21); TK Gristmill Software: Overland Park, KS, USA, 2014. [Google Scholar]
  41. Frisch, M.J.; Trucks, G.W.; Schlegel, H.B.; Scuseria, G.E.; Robb, M.A.; Cheeseman, J.R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G.A.; et al. Gaussian 09 D.01; Gaussian, Inc.: Wallingford, CT, USA, 2013. [Google Scholar]
  42. Verstraelen, T.; Pauwels, E.; de Proft, F.; van Speybroeck, V.; Geerlings, P.; Waroquier, M. Assessment of atomic charge models for gas-phase computations on polypeptides. J. Chem. Theory Comput. 2012, 8, 661–676. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  43. Verstraelen, T.; van Speybroeck, V.; Waroquier, M. The electronegativity equalization method and the split charge equilibration applied to organic systems: Parametrization, validation, and comparison. J. Chem. Phys. 2009, 131, 044127. [Google Scholar] [CrossRef] [PubMed]
  44. Seidler, T.; Champagne, B. Which charge definition for describing the crystal polarizing field and the χ(1) and χ(2) of organic crystals? Phys. Chem. Chem. Phys. 2015, 17, 19546–19556. [Google Scholar] [CrossRef] [PubMed]
  45. Martin, F.; Zipse, H. Charge distribution in the water molecule—A comparison of methods. J. Comput. Chem. 2005, 26, 97–105. [Google Scholar] [CrossRef] [PubMed]
  46. Rigby, J.; Izgorodina, E.I. Assessment of atomic partial charge schemes for polarisation and charge transfer effects in ionic liquids. Phys. Chem. Chem. Phys. 2013, 15, 1632–1646. [Google Scholar] [CrossRef] [PubMed]
  47. Bayly, C.I.; Cieplak, P.; Cornell, W.D.; Kollman, P.A. A well-behaved electrostatic potential based method using charge restraints for deriving atomic charges: The RESP model. J. Phys. Chem. 1993, 97, 10269–10280. [Google Scholar] [CrossRef]
  48. Wiberg, K.B.; Hadad, C.M.; LePage, T.J.; Breneman, C.M.; Frisch, M.J. Analysis of the effect of electron correlation on charge density distributions. J. Phys. Chem. 1992, 96, 671–679. [Google Scholar] [CrossRef]
  49. Stephens, M.E.; Becker, P.J. Virial partitioning analysis of electron correlation and nuclear motion in diatomic molecules. Mol. Phys. 1983, 49, 65–89. [Google Scholar] [CrossRef]
  50. Jabłoński, M.; Palusiak, M. Basis set and method dependence in atoms in molecules calculations. J. Phys. Chem. A 2010, 114, 2240–2244. [Google Scholar] [CrossRef] [PubMed]
  51. Jabłoński, M.; Palusiak, M. Basis set and method dependence in quantum theory of atoms in molecules calculations for covalent bonds. J. Phys. Chem. A 2010, 114, 12498–12505. [Google Scholar] [CrossRef] [PubMed]
  52. Davidson, E.R.; Chakravorty, S. A test of the Hirshfeld definition of atomic charges and moments. Theor. Chim. Acta 1992, 83, 319–330. [Google Scholar] [CrossRef]

Share and Cite

MDPI and ACS Style

Matczak, P. A Test of Various Partial Atomic Charge Models for Computations on Diheteroaryl Ketones and Thioketones. Computation 2016, 4, 3. https://doi.org/10.3390/computation4010003

AMA Style

Matczak P. A Test of Various Partial Atomic Charge Models for Computations on Diheteroaryl Ketones and Thioketones. Computation. 2016; 4(1):3. https://doi.org/10.3390/computation4010003

Chicago/Turabian Style

Matczak, Piotr. 2016. "A Test of Various Partial Atomic Charge Models for Computations on Diheteroaryl Ketones and Thioketones" Computation 4, no. 1: 3. https://doi.org/10.3390/computation4010003

APA Style

Matczak, P. (2016). A Test of Various Partial Atomic Charge Models for Computations on Diheteroaryl Ketones and Thioketones. Computation, 4(1), 3. https://doi.org/10.3390/computation4010003

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop