Approximate Reconstruction of Torsional Potential Energy Surface based on Voronoi Tessellation
1 Approximate Reconstruction of Torsional Potential Energy Surface based on Voronoi Tessellation
Chengming He, Yicheng Chi, and Peng Zhang*
Department of Mechanical Engineering The Hong Kong Polytechnic University, Hung Hom, Kowloon, Hong Kong
Corresponding Author:
Peng Zhang
Department of Mechanical Engineering
The Hong Kong Polytechnic University
Hung Hom, Kowloon, Hong Kong
E-mail: [email protected]
Fax: (852)23654703 Tel: (852)27666664
Colloquium:
Gas-phase Reaction Kinetics
Paper Length (Method 1):
The total word (excluding title block, abstract and the separate list of figure captions) is: words
Word Count (was performed from automatic counting function in MS Word plus Figures/Equations/References)
Abstract: 183 words, not included in word count
Main text : words Equations: 350 words
References: words (27 references)
Figures: words (8 figures with captions, color figures in electronic version only) Figure Column Word Count 1 Single 292 2 Single 172 3 Single 107 4 Double 346 5 Single 248 6 Single 118 7 Double 339 8 Single 241 Total words
One Supplemental Material is available. Approximate Reconstruction of Torsional Potential Energy Surface based on Voronoi Tessellation
Chengming He, Yicheng Chi, and Peng Zhang * Department of Mechanical Engineering The Hong Kong Polytechnic University, Hung Hom, Kowloon, Hong Kong
Abstract:
Torsional modes within a complex molecule containing various functional groups are often strongly coupled so that the harmonic approximation and one-dimensional torsional treatment are inaccurate to evaluate their partition functions. A family of multi-structural approximation methods have been proposed and applied in recent years to deal with the torsional anharmonicity. However, these methods approximate the exact “almost periodic” potential energy as a summation of local periodic functions with symmetric barrier positions and heights. In the present theoretical study, we illustrated that the approximation is inaccurate when torsional modes present non-uniformly distributed local minima. Thereby, we proposed an improved method to reconstruct approximate potential to replace the periodic potential by using information of the local minima and their Voronoi tessellation. First, we established asymmetric barrier heights by introducing two periodicity parameters and assuming that the exact barrier positions are at the boundaries of Voronoi cells. Second, we used multiplicatively weighted Voronoi tessellation to refine the barrier heights and positions by defining a structure-related distance metric. The proposed method has been tested for a few higher-dimensional cases, all of which show promising improved accuracy.
Keywords:
Partition function; Torsional anharmonicity; Multi-structural approximation; Voronoi tessellation; Almost periodic function; * Corresponding author E-mail address: [email protected] (P. Zhang) Fax: (852)23654703 Tel: (852)27666664 Accurate evaluation of partition functions is critical to thermochemical and kinetic calculations of complex molecules[1-4]. This is because torsional modes within a complex molecule containing various functional groups are often strongly coupled (SC), and the resulting torsional anharmonicity[5, 6] renders the harmonic approximation[7, 8] to be highly inaccurate. A variety of studies[9, 10] employ one-dimensional (1D) internal rotation treatment to replace a harmonic oscillator by a specific torsional mode, but it is often difficult to identify SC torsions with specific normal modes[11]. Then, it spawned the non-separable treatments of mixed torsions. The widely discussed Pitzer-Gwinn approximations[12], Feynman path integrals[13], and Monte Carlo phase space integrals[14] have been used but are usually computational expensive for complex molecules. In recent years, Zheng et al. [6, 15] proposed a family of multi-structural (MS) approximation methods, which were believed to satisfactorily deal with the torsional anharmonicity[16-19]. Among the MS methods, the MS-AS method[6] denoting “multi-structural method including all structures” is of particular interest because it does not require any information about conformational barriers or the paths that connect various structures. Furthermore, the exhaustive conformational structure search for all distinguishable structures is extremely complex and expensive for molecules with a large number of torsions (and hence a large number of conformational structures). Consequently, some cost-effective approximation methods, such as the MS-RS[6] method based on including all conformers generated from a reference structure (RS), a dual-level (a low-level and high-level electronic structure) method[20], and an extended two-dimensional (2D) torsion method[21] were also proposed. The MS-AS method[6] has been used in many studies[5, 17-19, 22]. For example, Zheng et al. [5] studied the hydrogen abstraction reactions of iso-butanol by hydroxyl radical, involving 9, 20, 18, 96, and 16 conformational structures for iso-butanol and the R1a−R1d transition states, respectively. They also showed the significance of considering the anharmonicity of high-frequency modes. Zheng et al. [22] studied the C–H bond dissociation processes of n-hexane and iso-hexane with containing 23 and 13 conformational structures in the parent molecules and 14 to 45 conformational structures in each of the seven isomeric products. Monge-Palacios et al. [17] studied the isomerization reaction of a six carbon atom Criegee intermediate (C6-CI) catalyzed by formic acid and found that the contribution of the reactant C6-CI conformers to the multi-structural partition function is larger than that of the saddle point conformers. Felsmaan et al. [19] applied the MS-AS method to predict partition functions of both reactants and transition states in MP (methyl propanoate) + OH/HO reaction systems, and they could successfully predict the rate constants compared with literature data. Shang et al. [18] used the MS-AS method to investigate the reaction kinetics of H-abstractions from DMA by H, CH , OH, and HO radicals in a broad temperature range (100–2000 K). We noted that the torsional anharmonicity results in an “almost periodic function” of potential energy on the parameter space of dihedral angles corresponding to the torsional modes. The MS-AS method invokes an essential approximation by expressing such a potential energy function as a summation of several local periodic functions. It means that the constructed potential energy surface (PES) along a specific dihedra angle is always symmetric on two sides of each local minimum. This approximation becomes increasingly inaccurate when torsional modes present non-uniformly distributed local energy minima. In the present study, we attempt to theoretically analyze the uncertainties of the MS-AS method and to propose improved methods for reconstructing approximate PES based on the mathematical technique of Voronoi tessellation. To deal with the torsional anharmonicity induced by SC torsional modes, the MS-AS method[6] uses Voronoi tessellation[6, 23] to identify the influence region for each conformational structure on the parameter space of dihedral angles so as to circumvent the difficulties of assigning coupled torsions to specific normal modes. As shown in Fig. 1, Voronoi tessellation cuts the entire parameter space into several one-by-one subspaces to each local minimum (solid points) and guarantees every point in a subspace is nearest to its corresponding local minimum than other minima. It can be mathematically described as that, for each structure (point) 𝑗 , its influence region 𝑅(𝑗) is defined as
𝑅(𝑗) = {𝑥 ∈ 𝑋|𝑑 (𝑥, 𝑗) ≤ 𝑑 (𝑥, 𝑘), for all 𝑗 ≠ 𝑘} (1) where 𝑥 is a point in the parameter space 𝑋 , and 𝑑 the distance metric. The Euclidean norm is generally used as the default distance metric. The periodic replicas are also included in the tessellation calculation so that could properly handle the periodic nature of the torsional modes. To physically reflect the influence region for each structure, a multiplicatively weighted Voronoi diagram[23] will be used in the present study. More details about the mathematical theory of Voronoi tessellation can be found in[23, 24]. Figure 1. Schematic (a) 2D and (b) 3D Voronoi tessellation generated from the presently calculated 78 conformational structures of MB (Methyl Butanoate) with HO in the transition state. Voronoi tessellation and reconstruction of PES require information from electronic structure calculations. In the present study, all conformational structures with local potential energy minima can be identified by internal torsions from an initial structure (generally a lowest energy state). Density functional theory (DFT) employing the B3LYP functional with the 6-311++G(d,p) basis set[25] was used for geometry optimization, frequency calculation, zero-point energy correction, and hindrance potential treatment. All the calculations were performed with the Gaussian 09 program package. 𝑄 CM = ( 12𝜋𝛽ℏ ) 𝑡 2⁄ (det{𝐃}) ∫ … ∫ 𝑑𝜙 … 𝑑𝜙 𝑡 𝑒 −𝛽𝑉(𝜙 ,…,𝜙 𝑡 )2𝜋 𝜎 𝑡 ⁄02𝜋 𝜎 ⁄0 (2) where 𝛽 = 1 𝑘 B 𝑇⁄ , 𝑘 B the Boltzmann’s constant, 𝑇 the temperature, ℏ the Planck’s constant, 𝐃 the torsional kinetic energy matrix that evaluated at the global minimum[27], 𝜙 𝜏 the torsional internal coordinate, 𝜎 𝜏 the parameter characterizing the periodicity of torsional space, and 𝑡 the total number of coupled torsions. Owing to the existence of a large number of local potential energy minima that contribute to the evaluation of partition function, the entire torsional space can be topographically divided into a number of distinct subspaces corresponding to each local minimum (structure), and thereby the total partition function is a summation of that for all the minima (structures). For each structure 𝑗 belonging to a certain torsion 𝜏 , the PES is assumed to be a periodic function: 𝑉 𝑗,𝜏 = 𝑈 𝑗 + 𝑊 𝑗,𝜏 𝑗,𝜏 (𝜙 𝜏 − 𝜙 𝜏.eq.𝑗 )], −𝜋𝑀 𝑗,𝜏 ≤ 𝜙 𝜏 − 𝜙 𝜏.eq.𝑗 ≤ 𝜋𝑀 𝑗,𝜏 (3) where 𝑈 𝑗 and 𝜙 𝜏,eq,𝑗 are respectively the potential energy and the torsional internal coordinate of the local minimum 𝑗 , 𝑀 𝑗,𝜏 the periodicity parameter. 𝑊 𝑗,𝜏 is the effective barrier height estimated from 𝑀 𝑗,𝜏 , the frequency 𝜔 𝑗,𝜏 , and the internal moment of inertia 𝐼 𝑗,𝜏 by 𝑊 𝑗,𝜏 = 2𝐼 𝑗,𝜏 𝜔 𝑗,𝜏2 𝑀 𝑗,𝜏2 (4) Furthermore, for the torsional anharmonicity induced by SC torsions, the PES within a certain structure 𝑗 is assumed to be separable so that 𝑉 𝑗 (𝜙 , … , 𝜙 𝑡 ) ≈ ∑ 𝑉 𝑗,𝜏 (𝜙 𝜏 ) 𝑡𝜏=1 (5) where each 𝑉 𝑗,𝜏 (𝜙 𝜏 ) is calculated by Eq. (3). The assignment of 𝑀 𝑗,𝜏 is a crucial part of the MS-AS method. For nearly separable (NS) torsions with approximately evenly distributed local minima, 𝑀 𝑗,𝜏 simply equals to the total number of local minima in the specific torsion; whereas for SC torsions, as shown in Fig. 1, 𝑀 𝑗,𝜏 is determined by Voronoi tessellation[6], in which 𝑀 𝑗,𝜏 is replaced by 𝑀 𝑗SC and assumed to be equivalent in every SC torsion. Then, the local periodicity 𝑀 𝑗,𝜏 is defined as 𝑀 𝑗,𝜏 = 𝑀 𝑗SC = 2𝜋(Ω 𝑗SC ) SC ⁄ (6) where Ω 𝑗SC is the hypervolume of subspace 𝑗 and 𝑡 SC is the total number of SC torsions. Overall, 𝑀 𝑗,𝜏 plays three important roles in the method. First, it controls the local periodicity. Second, it determines the integral subspace for a specific structure. Third, it accounts for the evaluation of implicit barrier heights. The periodic potential assumption of Eq. (3) in the MS-AS method is conditionally accurate if the local minima are uniformly distributed with same energies and frequencies. As the 1D example shown in Fig. 2, the potential 𝑉 cm −1 ⁄ = 750 cos(2𝜙) + 750 is a periodic function with the periodicity of 180 degree. Then it can be divided into two identical periodic functions that are corresponding to two local minima, respectively, in which each periodic function traverses a complete period of 180 degree. The exact potential is generally not a periodic function if existing the torsional anharmonicity, as the 1D torsional potential[6] for H O shown in Fig. 2, that the potential 𝑉 cm −1 ⁄ = 830.7 + is an almost periodic function with a permanent periodicity of 360 degree because of the natural torsion. The MS-AS method splits the potential 𝑉 into two piecewise periodic functions 𝑉 that corresponding to two local minima. According to Eq. (3) and (4) with the periodicity parameter 𝑀 = 2 , the internal moment of inertia
𝐼 = 0.4232 amu Å , and the harmonic frequency 𝜔 = 382.6 cm −1 for both two minima, the reconstructed potential 𝑉 is also plotted in Fig. 2. It is seen that the potential curve 𝑉 is always symmetric on two sides of the local minimum with same barrier heights, which causes the potential to be one side higher and the other side lower than the exact potential. It is also found an overlapped region between two neighboring periodic functions. Consequently, the periodic potential assumption properly leads to inaccurate estimation of the partition function, which is believed to be increasingly exacerbated when the local minima become more non-uniformly distributed. Figure 2. Potentials of a periodic function 𝑉 , an almost periodic function 𝑉 (corresponding to H O ), and piecewise periodic functions 𝑉 that constructed in the MS-AS method. The uncertainty of the MS-AS method becomes prominent when considering SC torsions. As the red frame of representative structures 1 and 2 shown in Fig. 1(a), the constructed subspace in the MS-AS method is a square for 2D and a cube for 3D coupled torsions. From a closeup of structure 1 shown in Fig. 3, we can infer that the Voronoi subspace and the constructed MS-AS subspace contribute approximately equally on the partition function for those structures with local minimum located almost centrally. However, for those structures whose local minima are very close to the boundary of their Voronoi cells, for example structure 2, the periodic potential assumption in the MS-AS method would induce highly inaccurate partition function estimation. This can be understood by a closeup of structure 2 shown in Fig. 3. Apart from the overlapped white region, the shade region denoted by vertical lines of isotropic subspace is closer to the local minimum than that denoted by horizontal lines of Voronoi subspace, and thereby has a lower potential energy 𝑉 and a high partition function 𝑄~𝑒 −𝑉 𝑘 B 𝑇⁄ . It is predicted that the deviations of partition function estimation between the MS-AS method and the exact value would be increasingly enlarged when considering three SC torsions in the transition state involving extremely non-uniformly distributed local minima. For example, Li et al. [16] found that the multi-structural torsional partition functions are 5.67 to 11.81 times larger than the single-structural partition functions for the transition state of the hydrogen abstraction reaction from methyl butenoate by H atoms. Figure 3. Closeup of Voronoi subspace and constructed MS-AS subspace in Fig. 1. 3.2 Improved reconstruction of torsional PES Considering the uncertainties of the MS-AS method, we considered to relax the assumption of the periodic potential and to obtain a physically reliable potential. We noted that a higher-level MS method (MS-ASCB)[6] has been introduced and ASCB denotes “based on all structures and conformational barriers”. It includes the explicit conformational barrier heights and positions obtained from electronic structure calculations, and hence the continuous torsional potential can be given as 𝑉 𝑗,𝜏 = { 𝑈 𝑗 + 𝑊 𝑗,𝜏𝐿 𝜏 − 𝜙 𝜏.eq.𝑗 )𝜋(𝜙 𝜏.eq.𝑗 − 𝜙 𝑗,𝜏𝐿 ))] , (𝜙 𝑗,𝜏𝐿 ≤ 𝜙 𝜏 ≤ 𝜙 𝜏.eq.𝑗 )𝑈 𝑗 + 𝑊 𝑗,𝜏𝑅 𝜏 − 𝜙 𝜏.eq.𝑗 )𝜋(𝜙 𝑗,𝜏𝑅 − 𝜙 𝜏.eq.𝑗 ))] , (𝜙 𝜏.eq.𝑗 ≤ 𝜙 𝜏 ≤ 𝜙 𝑗,𝜏𝑅 ) (7) where 𝑊 𝑗,𝜏𝐿 and 𝑊 𝑗,𝜏𝑅 are the exact left and right barrier heights along torsion 𝜏 of structure 𝑗 , and 𝜙 𝑗,𝜏𝐿 and 𝜙 𝑗,𝜏𝑅 are their corresponding barrier positions. Theoretically, the MS-ASCB method can provide more reliable potential than the MS-AS method. But, the additionally required information of barrier heights and positions would cause a significant amount of computational cost and human efforts. Consequently, it is quite compelling to think of a more physical potential to replace the periodic potential in the MS-AS method because the information of barrier positions and heights are reflected by the shape and boundary of the Voronoi tessellation. We attempt to reconstruct an almost periodic function of potential based on the Voronoi tessellation in two steps. First, we assume that the barrier positions are at the boundaries of Voronoi cells, and then calculate the barrier heights according to Eq. (4) by defining another two periodicity parameters, 𝑀 𝑗,𝜏𝐿 and 𝑀 𝑗,𝜏𝑅 , which are related to 𝑀 𝑗,𝜏 . As a result, we can characterize the asymmetric barrier positions and heights for the local minimum. Second, recognizing that barrier positions tend to vary with the energy and frequency of the local minimum, we define a structure-related distance metric to recalculate the Voronoi tessellation and to obtain the corrected barrier positions and heights. Specifically, the left and right boundaries of Voronoi cells in each torsion, 𝜙 𝑗,𝜏𝐿 and 𝜙 𝑗,𝜏𝑅 , are assumed as the exact barrier positions, which are already determined from the Voronoi tessellation, as the hollow intersection points of blue Voronoi subspace and red MS-AS subspace shown in Fig. 4(a). The local periodicity parameter 𝑀 𝑗,𝜏 for each structure is equivalent in every torsion and determined by Eq. (6). Then, to characterize the barrier height difference on the left and right sides of the local minimum in each torsion, we define two periodicity parameters, 𝑀 𝑗,𝜏𝐿 and 𝑀 𝑗,𝜏𝑅 , by 𝑀 𝑗,𝜏𝐿 Ω 𝑗,𝜏𝐿 = 𝑀 𝑗,𝜏𝑅 Ω 𝑗,𝜏𝑅 (8) (𝑀 𝑗,𝜏𝐿 ) + (𝑀 𝑗,𝜏𝑅 ) = 2𝑀 𝑗,𝜏2 (9) where Ω 𝑗,𝜏𝐿 and Ω 𝑗,𝜏𝑅 are the hypervolumes obtained by cutting the subspace into “left” and “right” parts with a specific plane, in which the specific plane is perpendicular to the torsion axis and across the internal coordinate 𝜙 𝑗,𝜏 of the local minimum, as shown in Fig. 4(a). Subsequently, we can obtain the corrected barrier heights, 𝑊 𝑗,𝜏𝐿 and 𝑊 𝑗,𝜏𝑅 , by Eq. (4), and the almost periodic potential by Eq. (7). Figure 4. A closeup of upper left corner of Fig. 1(a) showing the improvements of the MS-AS method with corrected (a) barrier heights based on original Voronoi diagram and (b) barrier positions based on multiplicatively weighted Voronoi diagram. In general cases, the barrier positions are not exactly at the boundaries of Voronoi cells, since the Voronoi tessellation topologically divides the space merely in geometry by assuming all structures contribute equally to the total potential. But in fact, the barrier positions would change as varying the energy and the frequency of each local minimum. Physically, for the ideal states of evenly distributed local minima with same energies and frequencies, the exact barrier positions should be completely located at the boundaries of Voronoi cells. However, for the general states involving frequency and energy differences between neighboring structures, the barrier positions tend to be closer to those minima with higher frequencies or higher energies. Based on this understanding, we define a new structure-related distance metric as 𝑑 𝑗,𝜏1 = 𝜀 𝑗,𝜏 𝑑 𝑗,𝜏0 , where 𝑑 𝑗,𝜏0 is the distance metric that used in the original Voronoi tessellation. The correction coefficient 𝜀 𝑗,𝜏 is defined as 𝜀 𝑗,𝜏 = 𝑑 𝑗,𝜏1 𝑑 𝜏0 = 𝑈 𝑗 + 2𝐼 𝑗,𝜏 𝜔 𝑗,𝜏2 𝑀 𝑗,𝜏2
1𝑁 ∑ (𝑈 𝑗 + 2𝐼 𝑗,𝜏 𝜔 𝑗,𝜏2 𝑀 𝑗,𝜏2 ) 𝑁𝑗=1 (10) to reflect the synergetic effects of frequency and energy changing on the variation of barrier positions. A weighted correction coefficient 𝜀 𝑗 is defined as the root mean square of all concerned torsions 𝜀 𝑗,𝜏 . Equation (10) shows that either increasing 𝑈 𝑗 or 𝜔 𝑗,𝜏 would result in an increased 𝜀 𝑗,𝜏 and 𝑑 𝑗,𝜏1 , and thereby a reduced hypervolume of Voronoi subspace 𝑗 . This can be clearly shown by the multiplicatively weighted Voronoi diagram in Fig. 4(b) that the Voronoi cells for structure 1’ ( 𝜀 = 0.5 ) and 2’ ( 𝜀 = 1.3 ) are enlarged and reduced, respectively. Furthermore, for any two neighboring local minima with same 𝑈 𝑗 and 𝜔 𝑗,𝜏 , their 𝜀 𝑗,𝜏 are the same (but not has to be 1) so that the original boundaries of Voronoi cells are not changed. We name the improved method as MS-ASB, where B denotes the corrected “barrier heights and positions” by information of Voronoi tessellation. It is noted that, for 1D case, the potential integration of Eq. (2) is performed in Voronoi subspace rather than the constructed isotropic subspace in the original MS-AS method. Whereas, for high dimensional cases, the potential integration is performed in a rectangle subspace, as the red frame shown in Fig. 4(a), rather than the Voronoi subspaces (arbitrary polygon) owing to the separability assumption of potential in each coupled torsion. To facilitate the application of the proposed MS-ASB method, the 1D example in a specific torsion is used to clearly present the derivation of asymmetric barrier heights and weighted Voronoi subspace in the Sec. 1 of the Supplemental Material, and then the high-dimensional coupled torsions can be thereby obtained in a straightforward manner. 3.3 Testing cases In the present study, the improved MS-ASB method was tested by a number of cases, such as 1D torsional potentials of H O and a toy model, ethanol radical, MB+HO , 1-pentyl radical, and 1-butanol radical. The process of partition function calculations and computational data used in the test cases have been included in the Sec. 2 and 3 of the Supplemental Material, respectively. The MATLAB code for these test cases are also attached so that interested readers can use it to reproduce all the results. Here, a 1D torsional potential of H O , a 2D torsional mode of 1-pentyl radical, and a 3D torsional mode of 1-butanol radical were chosen for better illustration. In the first testing case, the 1D torsion potential[6] of H O that has been discussed in Fig. 2 is used. Owing to the symmetric distribution of local minima, the exact barrier positions are just at the boundaries of Voronoi subspaces. As shown in Fig. 5(a), the reconstructed almost periodic function of potential by the MS-ASB method is more physically reliable than the symmetric potential by the MS-AS method. Figure 5(b) shows the partition function ratio (PFR) normalized by the exact value of partition function. The MS-AS method overestimates the partition function. However, the improved MS-ASB method can reduce the overestimation by a factor of 1.2 and obtain an accurate estimation of partition function that closer to the exact value particularly at high temperatures, regardless of some deviations at low temperature. This might be attributed to the slightly larger barrier height estimated by MS-ASB. Figure 5. Comparison of (a) the potential curve and (b) partition function ratio (PFR) of the 1D torsion potential[6] of H O . In the second testing case, the 2D torsional mode of 1-pentyl radical[6] with 15 conformational structures is used. Figure 6 shows the comparison of partition function Q
MS-ASB /Q MS-AS due to the lack of exact value. The calculation results by using the MS-ASCB method is however provided neither in[6] nor subsequent studies[15, 20, 21]. The ratio Q
MS-ASB /Q MS-AS is smaller than the unity, which is consistent with the theoretical analysis that the MS-AS method tend to predict a larger value of PF if there are many non-uniformly distributed local minima. This is because, geometrically, the weighted Voronoi tessellation in the MS-ASB method can modify the original “banded” Voronoi subspace in the MS-AS method to a more “round” subspace with the local minima located centrally, as shown in Fig. 7. Figure 6. Comparison of partition function Q
MS-ASB /Q MS-AS of the 2D torsional mode of 1-pentyl radical. Figure 7. Schematic of the (a) original Voronoi diagram in the MS-AS method and the (b) weighted Voronoi diagram in the MS-ASB method for the 2D torsional mode of 1-pentyl radical. In the third testing case, the 3D torsional mode of 1-butanol radical[6] with 29 conformational structures is used. It is seen that, as show in Fig. 8(a), the ratio Q
MS-ASB /Q MS-AS is close to the unity although there are more coupled torsions than that in Fig. 6. This is probably owing to the relatively uniform distribution of the local energy minima shown in Fig. 8(b). It implies that the improvement of the MS-ASB method mainly relies on the largely non-uniform distribution of the local minima as discussed in Fig. 7 rather than by considering more coupled torsions. Figure 8. Test of the 3D torsional mode of 1-butanol radical. (a) comparison of partition function Q
MS-ASB /Q MS-AS and (b) schematic of weighted Voronoi diagram. In summary, both the theoretical analysis and numerical results show that the MS-AS method tends to predict a larger value of PF, but the proposed MS-ASB method can reduce the possible overestimation by a factor that depends on the specific chemical systems. It is also worthy to mention that, for the 2D torsional mode of ethanol radical in the Supplemental Material, the ratio Q
MS-ASB /Q MS-AS is nearly the unity because the conformational structures are highly symmetric, which indicates that the MS-ASB method can be degenerated to the MS-AS method for the evenly distributed local minima so that to further consolidate the validations. Furthermore, the proposed MS-ASB method does not cause much additional cost because all the parameters used in it are obtained from the Voronoi Tessellation, which is the required procedure for other MS-based methods. These test cases in the present study show that the computational cost of the MS-ASB method is about twice of the MS-AS method by using our MATLAB code, which has been provided in the Supplemental Material. The most time-consuming process involved in the MS-ASB method is the calculation of 𝑀 𝑗,𝜏𝐿 and 𝑀 𝑗,𝜏𝑅 by the Monte Carlo sampling method. Motivated by the observations that the widely-used MS-AS method to deal with the torsional anharmonicity often overpredicts partition functions of complex molecular system, we performed an uncertainty analysis to show that its periodic potential assumption would result in inaccurate estimation of partition function when torsional modes present non-uniformly distributed local energy minima. To remedy this problem, we proposed the improved MS-ASB method to reconstruct an almost periodic function of potential to replace the periodic potential based on the Voronoi tessellation of local minima. The testing cases in the present study show that the MS-ASB method is promising to predict more accurate PF compared with the MS-AS method when addressing the torsional anharmonicity problem with extremely non-uniformly distributed local minima. The proposed MS-ASB method can be treated as an intermediate method, which aims to have the comparable computational cost to the MS-AS method but to have computational accuracy approximate to that of the MS-ASCB method. Separate future works are merited to further validate the MS-ASB method against the MS-ASCB method and other higher-level theoretical methods, which often require the very large amount of computational costs for exact barrier heights and positions of couple torsional modes.
Acknowledgment This work was supported by NSFC (91641105) and partly by the university matching grant (4-BCE8) and DGRF (G-UAHP).
References [1]. Klippenstein, S.J., V.S. Pande, and D.G. Truhlar, Chemical kinetics and mechanisms of complex systems: a perspective on recent theoretical advances. J. Am. Chem. Soc., 2014. (2): p. 528-546. [2]. Strekalov, M., An accurate closed-form expression for the partition function of Morse oscillators. Chem. Phys. Lett., 2007. (1-3): p. 209-212. [3]. Beste, A., One-dimensional anharmonic oscillator: Quantum versus classical vibrational partition functions. Chem. Phys. Lett., 2010. (1-3): p. 200-205. [4]. Jia, C.-S., et al., Partition function of improved Tietz oscillators. Chem. Phys. Lett., 2017. : p. 150-153. [5]. Zheng, J., R.n. Meana-Pañeda, and D.G. Truhlar, Prediction of experimentally unavailable product branching ratios for biofuel combustion: The role of anharmonicity in the reaction of isobutanol with OH. J. Am. Chem. Soc., 2014. (13): p. 5150-5160. [6]. Zheng, J., et al., Practical methods for including torsional anharmonicity in thermochemical calculations on complex molecules: The internal-coordinate multi-structural approximation. Phys. Chem. Chem. Phys, 2011. (23): p. 10885-10907. [7]. Ayala, P.Y. and H.B. Schlegel, Identification and treatment of internal rotation in normal mode vibrational analysis. J. Chem. Phys., 1998. (6): p. 2314-2325. [8]. Lin, C.Y., E.I. Izgorodina, and M.L. Coote, How accurate are approximate methods for evaluating partition functions for hindered internal rotations? J. Phys. Chem. A, 2008. (9): p. 1956-1964. [9]. Ellingson, B.A., et al., Statistical thermodynamics of bond torsional modes: Tests of separable, almost-separable, and improved Pitzer–Gwinn approximations. J. Chem. Phys., 2006. (8): p. 084305. [10]. Sharma, S., S. Raman, and W.H. Green, Intramolecular hydrogen migration in alkylperoxy and hydroperoxyalkylperoxy radicals: accurate treatment of hindered rotors. J. Phys. Chem. A, 2010. (18): p. 5689-5701. [11]. Van Speybroeck, V., D. Van Neck, and M. Waroquier, Ab initio study of radical reactions: role of coupled internal rotations on the reaction kinetics (III). J. Phys. Chem. A, 2002. (38): p. 8945-8950. [12]. Pitzer, K.S. and W.D. Gwinn, Energy levels and thermodynamic functions for molecules with internal rotation: I. Rigid frame with attached tops, in Mol. Struct. Stat. Thermodyn. 1993, World Scientific. p. 33-46. [13]. Mielke, S.L. and D.G. Truhlar, Improved methods for Feynman path integral calculations and their application to calculate converged vibrational–rotational partition functions, free energies, enthalpies, entropies, and heat capacities for methane. J. Chem. Phys., 2015. (4): p. 044105. [14]. Jasper, A.W., et al., Anharmonic Rovibrational Partition Functions for Fluxional Species at High Temperatures via Monte Carlo Phase Space Integrals. J. Phys. Chem. A, 2018. (6): p. 1727-1740. [15]. Zheng, J. and D.G. Truhlar, Quantum thermochemistry: Multistructural method with torsional anharmonicity based on a coupled torsional potential. J. Chem. Theory Comput., 2013. (3): p. 1356-1367. [16]. Li, X., et al., Kinetics and branching fractions of the hydrogen abstraction reaction from methyl butenoates by H atoms. Phys. Chem. Chem. Phys, 2017. (25): p. 16563-16575. [17]. Monge-Palacios, M., et al., Theoretical kinetic study of the formic acid catalyzed Criegee intermediate isomerization: multistructural anharmonicity and atmospheric implications. Phys. Chem. Chem. Phys, 2018. (16): p. 10806-10814. [18]. Shang, Y., et al., Chemical kinetics of H-abstractions from dimethyl amine by H, CH3, OH, and HO2 radicals with multistructural torsional anharmonicity. Phys. Chem. Chem. Phys, 2019. (23): p. 12685-12696. [19]. Felsmann, D., et al., Contributions to improving small ester combustion chemistry: Theory, model and experiments. Proc. Combust. Inst., 2017. (1): p. 543-551. [20]. Bao, J.L., L. Xing, and D.G. Truhlar, Dual-Level Method for Estimating Multistructural Partition Functions with Torsional Anharmonicity. J. Chem. Theory Comput., 2017. (6): p. 2511-2522. [21]. Simón-Carballido, L., et al., Anharmonicity of Coupled Torsions: The Extended Two-Dimensional Torsion Method and Its Use To Assess More Approximate Methods. J. Chem. Theory Comput., 2017. (8): p. 3478-3492. [22]. Zheng, J., T. Yu, and D.G. Truhlar, Multi-structural thermodynamics of C–H bond dissociation in hexane and isohexane yielding seven isomeric hexyl radicals. Phys. Chem. Chem. Phys, 2011. (43): p. 19318-19324. [23]. Mu, L., Polygon characterization with the multiplicatively weighted Voronoi diagram. Prof. Geogr., 2004. (2): p. 223-239. [24]. Aurenhammer, F., Voronoi diagrams—a survey of a fundamental geometric data structure. ACM Comput. Surv. , 1991. (3): p. 345-405. [25]. Becke, A.D., Density-functional thermochemistry. III. The role of exact exchange. J. Chem. Phys., 1993. (7): p. 5648. [26]. Kilpatrick, J.E. and K.S. Pitzer, Energy levels and thermodynamic functions for molecules with internal rotation. III. Compound rotation. J. Chem. Phys., 1949. (11): p. 1064-1075. [27]. Zhang, L. and P. Zhang, Towards high-level theoretical studies of large biodiesel molecules: an ONIOM [QCISD (T)/CBS: DFT] study of hydrogen abstraction reactions of C n H 2n+ 1 COOC m H 2m+ 1+ H. Phys. Chem. Chem. Phys, 2015. (1): p. 200-208. List of Figure captions (Color figures in electronic version only) Figure 1. Schematic (a) 2D and (b) 3D Voronoi tessellation generated from the presently calculated 78 conformational structures of MB (Methyl Butanoate) with HO in the transition state. Figure 2. Potentials of a periodic function 𝑉 , an almost periodic function 𝑉 (corresponding to H O ), and piecewise periodic functions 𝑉 that constructed in the MS-AS method. Figure 3. Closeup of Voronoi subspace and constructed MS-AS subspace in Fig. 1. Figure 4. A closeup of upper left corner of Fig. 1(a) showing the improvements of the MS-AS method with corrected (a) barrier heights based on original Voronoi diagram and (b) barrier positions based on multiplicatively weighted Voronoi diagram. Figure 5. Comparison of (a) the potential curve and (b) partition function ratio (PFR) of the 1D torsion potential[6] of H O Figure 6. Comparison of partition function Q
MS-ASB /Q MS-AS of the 2D torsional mode of 1-pentyl radical. Figure 7. Schematic of the (a) original Voronoi diagram in the MS-AS method and the (b) weighted Voronoi diagram in the MS-ASB method for the 2D torsional mode of 1-pentyl radical. Figure 8. Test of the 3D torsional mode of 1-butanol radical. (a) comparison of partition function Q
MS-ASB /Q MS-AS and (b) schematic of weighted Voronoi diagram. Supplemental Material
Item File name Content (captions) 1 PROCI-D-19-01430R2-Supplemental Materials.docx 1. Derivation of asymmetric barrier heights and weighted Voronoi subspace 錯 誤 ! 尚未定義書籤。
2. Process of partition function calculations by attached MATLAB codes .. 錯誤 ! 尚未定義書籤。 Run_1_preperation.m executable 錯誤 ! 尚未定義書籤。 Run_2_weighted_voronoi.m executable 錯誤 ! 尚未定義書籤。 Run_3_pf.m executable ............. 錯誤 ! 尚未定義書籤。 executable ............... 錯誤 ! 尚未定義書籤。
3. Computational data for testing cases ......................... 錯誤 ! 尚未定義書籤。 O ................. 錯誤 ! 尚未定義書籤。 錯誤 ! 尚未定義書籤。 錯誤 ! 尚未定義書籤。 radical 錯 誤 ! 尚 未 定 義 書籤。 錯誤 ! 尚未定義書籤。 錯 誤 ! 尚 未 定 義 書籤。
4. Summary ........................................................................ 錯誤 ! 尚未定義書籤。 Reference ............................................................................ 錯誤 ! 尚未定義書籤。
2 PROCI-D-19-01430R2-matlab_codes.zip It contains 7 folders: /matlab_codes/H2O2 /matlab_codes/artificial_model /matlab_codes/ethanol /matlab_codes/MB+HO2 /matlab_codes/1-pentyl-2d /matlab_codes/1-pentyl-3d /matlab_codes/1-butanol-3d Each folder has four executables: