- Email: [email protected]
Calculation of packing structure of methanol solid using ab initio lattice energy at the MP2 level
Chemical Physics Letters 369 (2003) 597–604 www.elsevier.com/locate/cplett
Calculation of packing structure of methanol solid using ab initio lattice energy at the MP2 level Kanade Nagayoshi a,b, Kazuo Kitaura c,*, Shiro Koseki a, Suyong Re b, Kaoru Kobayashi b, Yoong-Kee Choe b, Shigeru Nagase b a
c
College of Integrated Arts and Sciences, Osaka Prefecture University, Sakai, Osaka 599-8531, Japan b Department of Theoretical Studies, Institute for Molecular Science, Okazaki 444-8585, Japan Research Institute for Computational Sciences, AIST, 1-1-1 Umezono, Tsukuba, Ibaraki 305-8568, Japan Received 20 November 2002; in final form 9 December 2002
Abstract The ab initio MO based lattice energy minimization method, which have been recently proposed, is applied to the calculation of the packing structures of methanol crystal, a- and b-phases. We employed the correlated wavefunction level of theory, the second order Møller–Plesset perturbation theory (MP2), for the calculation of the lattice energy. The lattice parameters optimized at the MP2/6-31++G** level for a-phase were in good agreement with the experimental values and these for b-phase were in moderate agreement with the experiment. The averaged relative errors of the calculated cell lengths were 0.9% and 6.2% for a- and b-phases, respectively. Ó 2003 Elsevier Science B.V. All rights reserved.
1. Introduction One of the challenging issues in physical chemistry is to develop a method for predicting crystal structures from the information of only their chemical compositions [1]. The major approaches for the molecular crystal structure calculations are based on statistical mechanics; Monte Carlo and molecular dynamics simulations have been performed using empirical potential functions [2,3]. Despite their advantage in computational costs, there is much difficulty to gener-
*
Corresponding author. Fax: +81-298-61-3171. E-mail address: [email protected] (K. Kitaura).
ate reliable potential functions. The other approach is to obtain the crystal structure via the lattice energy minimization method which costs less than statistical methods. However, it is difficult to calculate the total energy of crystal with the ab initio MO method from the point of computational costs, although the performance of computer is progressing. We have proposed a computational procedure to obtain the lattice energy directly from ab initio MO calculations without using empirical potential functions such as Lennard–Jones and Buckingham potentials [4]. In this method, the difficulty in generating parameters of potential functions is replaced by simple selection of basis set and wavefunction. Therefore, a variety of molecular
0009-2614/03/$ - see front matter Ó 2003 Elsevier Science B.V. All rights reserved. doi:10.1016/S0009-2614(03)00025-3
598
K. Nagayoshi et al. / Chemical Physics Letters 369 (2003) 597–604
crystals can be handled with our method. The method was first applied to the packing structure calculation of electron donor–acceptor complex H3 N–BF3 at the HF level. In this system, the hydrogen bonding between the complexes plays a major role to determine the packing structure. For non-polar molecules and molecules with non-polar groups, it is necessary to take electron correlation effect into account in the calculation of lattice energy for a reasonable description of solid structures. In this Letter, to assess the validity of the method at the correlated wavefunction level, we calculated the packing structure of methanol molecule that has a non-polar CH3 group. For methanol solid, three phases have been known; ordered a-phase, disordered b-phase, and metastable c-phase [5–7]. The former two phases have hydrogen bonded zigzag chain structures. The structure of c-phase is in between a- and b-phases. For a-phase, the molecular simulations with potential functions have been performed [2,3,8–10]. These molecular simulation studies have suggested necessity for construction of potential functions based on ab initio potential surfaces.
2. Method The detailed description of the computational procedure of the ab initio MO based lattice energy calculation is given elsewhere [4]. Here, the outline of the method is described. The lattice energy, Elatt , is given as follows: ( ) X X1 eij ; Elatt ¼ ei þ ð1Þ 2 i j6¼i where i runs over all molecules within an Ôasymmetric unitÕ in which the molecule is regarded as an unit object for symmetry operations associated with the space group of crystal. j runs over all the other molecules generated by the symmetry operations. The one-body energy ei is the polarization energy of the molecule i in crystal. The two-body energy eij is the intermolecular interaction energy between molecules i and j. In our method, these energies are directly obtained from ab initio MO calculations.
In this work, the MP2 method is applied to lattice energy calculations. The polarization energy is still evaluated at the HF level and only the pair interaction energy is calculated at the MP2 level using the following equation: n o polðjÞ polðiÞ ei þ ej ; eij ¼ EijMP2 EiMP2 þ EjMP2 ð2Þ where E with the superscript MP2 denotes the total energy of molecule or molecular pair at the MP2 level. Then, the first term of Eq. (2) is the usual intermolecular interaction energy between molecules i and j obtained from supermolecule calcupolðjÞ lations. ei in Eq. (2) is the polarization energy of molecule i involved in the interaction energy between molecules i and j and is calculated at the HF level by the same manner as the one-body polarization energy [4]. 3. Computational details The lattice energy minimization calculations were performed on methanol crystal at the MP2 level with various basis sets; 3-21G, 3-21+G, 3-21G**, 6-31G, 6-31+G, 6-31G*, 6-31+G*, 6-31++G, 6-31G**, and 6-31++G**. We focus on a- and b-phases which have zigzag chains of hydrogen bonds along b-axis in the former and c-axis in the latter (Fig. 1). In the calculations, we assumed ordered structure for both a- and b-phases. Since c-phase is metastable and its structure is unknown [6], no consideration is given to this phase. In order to avoid a multi-minimum problem, the initial parameters including the lattice constants, the molecular orientation, space group, and the number of molecules in unit cell (Z) were taken from the experimental values [5,6] (vide infra). The lattice energy was then minimized with respect to the lattice constants (a, b, and c) and the molecular orientation angles using numerical differentiation. Each methanol molecule in the crystal was kept rigid at the gas phase structure obtained by ab initio MO calculations at the same level. In the lattice energy calculations, we considered from asymmetric unit molemolecules within 20 A cules. The number of molecules within this cut-off length were 664 and 614 for a- and b-phases,
K. Nagayoshi et al. / Chemical Physics Letters 369 (2003) 597–604
599
Fig. 1. The structure of methanol crystals: (a) a-phase and (b) b-phase. The dashed lines show the unit cell.
respectively. Comparing with the lattice energy obtained from the calculation with larger cut-off , the number of methanol molecules: 10 434 for (50 A a-phase and 9606 for b-phase), we estimated that the error due to this approximation is 0.00 and 0.11 kcal/mol for a- and b-phases, respectively, at the MP2/6-31++G** level. Since the error in the lattice energy for a-phase was less than 0.01 kcal/ mol, the Ewald summation [11] was not used. To make even the calculation level and to make sure the difference of phases, the Ewald summation was purposely not used in the calculation of b-phase. The following is the summary of the initial parameters used for the lattice energy minimizations. For a-phase; orthorhombic, P 21 21 21 space group, , and Z ¼ 4 at a ¼ 4:873, b ¼ 4:461, c ¼ 8:867 A 15 K [6]. For b-phase; orthorhombic, Cmcm space , and Z ¼ 4 group, a ¼ 6:43, b ¼ 7:24, c ¼ 4:67 A at 165 K [5]. All calculation were performed with our original program which utilizes the GA U S S I A N 94 [12].
4. Results and discussion For a basis to understand the molecular interaction in the crystal, the optimized structures of
the two isomers of methanol dimer obtained from MP2/6-31++G** calculations are shown in Fig. 2. One is Ôhead to headÕ type (hh) and the other is Ôhead to tailÕ type (ht). The ÔheadÕ means the hydroxyl group and the ÔtailÕ means the methyl group. hh has the linear hydrogen bond between each hydroxyl group and ht has the cyclic hydrogen bond through the two OH groups and the methyl groups. In hh, the O HðOÞ hydrogen and the O–O distance is bond length is 1.888 A . The OH O part in hh is almost linear 2.854 A (\OH O ¼ 172:4°). On the other hand, in ht, there are two weak O HðCÞ type hydrogen bond. The O HðCÞ distance is significantly long ) and the \CH O angle (106.1) is lar(2.401 A gely deviated from 180°. Thus, the Ôhead to tailÕ type of the cyclic hydrogen bonding interaction between two methanol molecules is very weak, even though the hydrogen bond is multiple. For reason of those, the lowest energy configuration of the dimer is linear rather than cyclic. The results are in good agreement with the previous theoretical study [13]. The intermolecular interaction energy of hh ()7.2 kcal/mol) is more stable than that of ht ()3.3 kcal/mol) by 3.9 kcal/mol.This suggests that the Ôhead to headÕ type of hydrogen bonding interaction dominates the crystal packing.
600
K. Nagayoshi et al. / Chemical Physics Letters 369 (2003) 597–604
and Fig. 2. The optimized structures of ðCH3 OHÞ2 at the MP2/6-31++G** level of theory. The bond lengths and angles are in A degrees, respectively. The intermolecular interaction energies, E, are in kcal/mol.
and 11.5°, respectively. The avangles are 0.49 A eraged relative error in the cell lengths is 6.2%. These errors are non-negligibly larger than those of a-phase. This result of b-phase can not be compared with the force field calculations, since they have not been reported. The calculated lattice constants and the molecular orientation highly depend on the basis set employed. We examined the basis set dependence of these quantities by using various basis sets (Fig. 3). The poorest result was obtained at the MP2/3-21G* level; the averaged relative error in the cell lengths exceeded 10% for a-phase. Nevertheless, the order of cell lengths (c > a > b) is well reproduced at most of the basis sets examined. Similar result was also obtained for b-phase.
Table 1 shows the optimized lattice constants of a- and b-phases at the MP2/6-31++G** level. For a-phase, both the lattice constants and the molecular orientation are in good agreement with the experimental values. The maximum deviation of calculated values from the experimental ones [6] is for the lattice constants and 2.2° for the 0.08 A angles of the molecular orientation. The averaged relative error in cell lengths is 0.9%. Note that the results of the force field calculation reported by Mooij et al. [2] are in good agreement with the present results. For b-phase, the optimized lattice constants and the molecular orientation at the MP2/6-31++G** level are in moderate agreement with the experimental values [5]; the maximum deviations of the cell lengths and the orientation
Table 1 Calculated lattice constants and molecular orientation of methanol solid, a- and b-phases, at the MP2/6-31++G** levela a-phase
a b c w h u Lattice energy
b-phase b
c
This work
Force field
Obs
4.85 ()0.02) 4.56 ()0.08) 8.83 ()0.04) 65.4 (+0.2) 85.4 (+2.2) 26.7 (+0.4) )15.0
5.07 (+0.20) 4.66 (+0.02) 8.71 ()0.16) – – – –
4.873 4.641 8.867 65.2 83.2 26.3 –
This work
Obsd
6.08 ()0.35) 7.43 (+0.19) 5.16 (+0.49) 78.5 ()11.5) 185.6 (+5.6) 85.7 ()4.3) )11.4
6.43 7.24 4.67 90.0 180.0 90.0 –
and degrees, respectively. The data in The lattice constants (a; b, and c) and the molecular orientations (w; h, and u) are in A brackets are deduction from observed values. The angles w; h, and u are defined as follows: the rotation by w; h, and u about the molecule lying in the ac plane and C–O axis along the c-axis is rotated to the observed position by w; h, and u about the a-, b-, and c-axes, respectively. b Ref. [2]. The AMMs-II force filed. c Ref. [6]. d Ref. [5]. a
K. Nagayoshi et al. / Chemical Physics Letters 369 (2003) 597–604
Fig. 3. The averaged relative error of the calculated cell lengths. The error is given by fja0 aj=a þ jb0 bj=b þ jc0 cj=cg=3 where a; b, and c are the experimental values and a0 ; b0 , and c0 are calculated ones.
As in the case of a-phase, there is large basis set dependence of the lattice constants, but the order of cell lengths (b > a > c) is again well predicted at most of the basis sets. Although the basis set dependence appears to be not consistent, inclusion of diffuse functions generally reduces the averaged relative error in the cell lengths by 0.5% on average. It is noted that the rather small 321+G basis set gives a good results; the MP2/321+G calculations give the lattice constants in comparable accuracy to the MP2/6-31++G** calculations. Now, let us look at the pair interactions of the molecules in crystal together with their geometries obtained from the MP2/6-31++G** calculation. To simply indicate molecular pairs, they are named as follows; for example, the pair of molecules 1 and 2 (Fig. 1) with the Ôhead to headÕ type interaction in a-phase is called A12-hh, that in b-zphase B12-hh, and so forth. Fig. 4a shows the important molecular pairs in a-phase, whose stabilization energy is more than 1 kcal/mol. Clearly, the most important molecular pair is A12-hh which corresponds to the structure of the most stable dimer (hh) in gas phase. Due to the packing effect, the O–O distance of A12-hh is significantly short
601
) compared to that of the dimer hh in gas (2.697 A ). The other molecular pairs, A13phase (2.854 A ht, A14-ht, and A15-ht, gain rather small stabilization through the distorted Ôhead to tailÕ type interaction. The representative molecular pairs in b-phase are shown in Fig. 4b. A dominant molecular pair is B12-hh which has the Ôhead to headÕ type hydrogen bonding. This pair is resemble to that of A12-hh. As in the case of a-phase, again, one could see the significant shortening of the O–O ) owing to the packing effect, distance (2.740 A is smaller than though the decrease of 0.114 A in a-phase. On the other hand, that of 0.157 A the stabilization energy of the molecular pair B13ht is small due to the weak interaction of Ôhead to tailÕ type. A12-hh and B12-hh form hydrogen bond chain along a certain axis in the crystal. In a-phase, the chain grows along the b-axis (Fig. 5a). The calcu is in reasonable lated O–O distance of 2.697 A . agreement with the experimental value of 2.74 A In b-phase, hydrogen bond chain grows along c-axis (Fig. 5b). The calculated O–O distance of is slightly longer than that in a-phase. On 2.740 A the other hand, the X-ray diffraction study has revealed that the O–O distance of the hydrogen ) is much shorter bond chains in b-phase (2.66 A ). This fact indicates than that of a-phase (2.74 A that there is a significant difference in the packing effect between a- and b-phases. For b-phase, the chain interaction of the Ôhead to headÕ type of hydrogen bonding along c-axis plays a dominant role in the packing. For a-phase, the several weak Ôhead to tailÕ type interactions along a- and c-axes are also important as well as the chain interaction of hydrogen bond along b-axis. Table 2 shows the calculated lattice energy with its component energies and the pair interaction energy of some nearest neighbors. The lattice energy of a-phase ()15.0 kcal/mol) is larger in magnitude than that of b-phase ()11.4 kcal/mol). Although no lattice energy has been observed experimentally, this result is consistent with the fact that a-phase is stable at lower temperature than b-phase. The polarization energy of a-phase is 1.5 kcal/mol larger in magnitude than that of bphase. Similarly, the Coulomb interaction energy
602
K. Nagayoshi et al. / Chemical Physics Letters 369 (2003) 597–604
Fig. 4. The geometrical parameters of molecular pairs in the optimized crystal structure at the MP2/6-31++G** level. The bond and degrees, respectively. The intermolecular interaction energies, E, are in kcal/mol. (a) is a-phase and lengths and angles are in A (b) is b-phase.
of a-phase is 4.4 kcal/mol larger in magnitude than that of b-phase. These results clearly show that the stabilization of a molecule due to surroundings is larger in a-phase than in b-phase. In a-phase, the interaction energy of the dominant molecular pair (A12-hh) provides only 33% of the whole lattice energy. In other word, a large portion of the lattice energy is brought by the many weak interactions. On the other hand, in b-phase, more than 45% of the lattice energy is provided by the dominant
molecular pair (B12-hh). These results also support that there is an apparent difference in nature of the crystal packing between a- and b-phases.
5. Summary We calculated the packing structure of the hydrogen bonding methanol crystal using the ab initio MO based lattice energy minimization
K. Nagayoshi et al. / Chemical Physics Letters 369 (2003) 597–604
603
) in hydrogen bonded chain in the methanol crystal. In the bracket is given the experimental value. Fig. 5. The O–O distance (A (a) a-phase and (b) b-phase.
Table 2 Calculated lattice energy and the intermolecular interaction energies of some nearest neighborsa
Lattice energy Polarization energy Coulomb energy Correlation energy Pair energyb 1-2 1-3 1-4 1-5 1-6 1-7 1-8
a-phase
b-phase
)15.0 )3.5 )14.6 )9.0
)11.4 )2.0 )10.2 )6.7
)4.9 )1.2 )1.1 )1.0 )0.7 )0.5 )0.4
(2) (2) (2) (2) (2) (2) (2)
)5.1 )1.0 )1.0 )0.7 )0.6 )0.2 0.1
(2) (2) (2) (2) (2) (2) (2)
constants, the molecular orientation, and hydrogen bond network in the crystal with satisfactory accuracy. The calculation also revealed the difference in the packing feature between a- and b-phases. Since the examined two phases have different nature of interaction, the simple point charge approximation employed in the estimation of the polarization energy may result in an unbalanced error in the lattice constants. We believe that our computational procedure for the lattice energy would be improved by using the fragment MO method [14,15] which gives the accurate total energy of molecular clusters based on a pair interaction approximation.
a
All energies are given in kcal/mol. Pair energies calculated by Eq. (2). See Fig. 1 for the indices of molecular pairs. The number of the equivalent pairs is given in parenthesis. b
method at the correlated wavefunction level. The calculated lattice constants depended highly on the basis set employed. It has been shown that the MP2/6-31++G** calculations predicts the lattice
Acknowledgements The numerical calculations were partially carried out at the Computer Center of the Institute for Molecular Science. This work was partly supported by ACT-JST.
604
K. Nagayoshi et al. / Chemical Physics Letters 369 (2003) 597–604
References [1] John Maddox, Nature 335 (1988) 201. [2] W.T.M. Mooij, B.P. van Eijck, J. Kroon, J. Phys. Chem. A. 103 (1999) 9883. [3] M. Hloucha, A.K. Sum, S.I. Sandler, J. Chem. Phys. 113 (2000) 5401. [4] T. Ikeda, K. Nagayoshi, K. Kitaura, Chem. Phys. Lett., in press. [5] K.J. Tauer;, W.N. Lipscomb, Acta Cryst. 5 (1952) 606. [6] B.H. Torrie, S.-X. Weng, B.M. powell, Mol. Phys. 67 (1989) 575. [7] B.H. Torrie, O.S. binbrek, M. Strauss, I.P. Swainson, J. Solid State Chem. 166 (2002) 415. [8] A.K. Sum, S.I. Sandler, R. Bukowski, K. Szalewicz, J. Chem. Phys. 116 (2002) 7627.
[9] R. Bukowski, K. Szalewicz, C.F. Chabalowski, J. Phys. Chem. A. 103 (1999) 7322. [10] J. Alonso, F.J. Bermejo, M. Garcia-Hernandez, J.L. Martines, W.S. Howells, A. Criado, J. Chem. Phys. 96 (1992) 7696. [11] P.P. Ewald, Ann. Physik. 64 (1921) 253. [12] M.J. Frisch et al., GA U S S I A N 94, Revision D.4, Gaussian, Pittsburgh, PA, 1995. [13] R.A. Provencal, J.B. Paul, K. Roth, C. Chapo, R.N. Casaes, R.J. Saykally, G.S. Tschumper, H.F. Schaefer III, J. Chem. Phys. 110 (1999) 4258. [14] K. Kitaura, T. Sawai, T. Asada, T. Nakano, M. Uebayashi, Chem. Phys. Lett. 312 (1999) 319. [15] K. Kitaura, E. Ikeo, T. Asada, T. Nakano, M. Uebayash, Chem. Phys. Lett. 313 (1999) 701.
Calculation of packing structure of methanol solid using ab initio lattice energy at the MP2 level Kanade Nagayoshi a,b, Kazuo Kitaura c,*, Shiro Koseki a, Suyong Re b, Kaoru Kobayashi b, Yoong-Kee Choe b, Shigeru Nagase b a
c
College of Integrated Arts and Sciences, Osaka Prefecture University, Sakai, Osaka 599-8531, Japan b Department of Theoretical Studies, Institute for Molecular Science, Okazaki 444-8585, Japan Research Institute for Computational Sciences, AIST, 1-1-1 Umezono, Tsukuba, Ibaraki 305-8568, Japan Received 20 November 2002; in final form 9 December 2002
Abstract The ab initio MO based lattice energy minimization method, which have been recently proposed, is applied to the calculation of the packing structures of methanol crystal, a- and b-phases. We employed the correlated wavefunction level of theory, the second order Møller–Plesset perturbation theory (MP2), for the calculation of the lattice energy. The lattice parameters optimized at the MP2/6-31++G** level for a-phase were in good agreement with the experimental values and these for b-phase were in moderate agreement with the experiment. The averaged relative errors of the calculated cell lengths were 0.9% and 6.2% for a- and b-phases, respectively. Ó 2003 Elsevier Science B.V. All rights reserved.
1. Introduction One of the challenging issues in physical chemistry is to develop a method for predicting crystal structures from the information of only their chemical compositions [1]. The major approaches for the molecular crystal structure calculations are based on statistical mechanics; Monte Carlo and molecular dynamics simulations have been performed using empirical potential functions [2,3]. Despite their advantage in computational costs, there is much difficulty to gener-
*
Corresponding author. Fax: +81-298-61-3171. E-mail address: [email protected] (K. Kitaura).
ate reliable potential functions. The other approach is to obtain the crystal structure via the lattice energy minimization method which costs less than statistical methods. However, it is difficult to calculate the total energy of crystal with the ab initio MO method from the point of computational costs, although the performance of computer is progressing. We have proposed a computational procedure to obtain the lattice energy directly from ab initio MO calculations without using empirical potential functions such as Lennard–Jones and Buckingham potentials [4]. In this method, the difficulty in generating parameters of potential functions is replaced by simple selection of basis set and wavefunction. Therefore, a variety of molecular
0009-2614/03/$ - see front matter Ó 2003 Elsevier Science B.V. All rights reserved. doi:10.1016/S0009-2614(03)00025-3
598
K. Nagayoshi et al. / Chemical Physics Letters 369 (2003) 597–604
crystals can be handled with our method. The method was first applied to the packing structure calculation of electron donor–acceptor complex H3 N–BF3 at the HF level. In this system, the hydrogen bonding between the complexes plays a major role to determine the packing structure. For non-polar molecules and molecules with non-polar groups, it is necessary to take electron correlation effect into account in the calculation of lattice energy for a reasonable description of solid structures. In this Letter, to assess the validity of the method at the correlated wavefunction level, we calculated the packing structure of methanol molecule that has a non-polar CH3 group. For methanol solid, three phases have been known; ordered a-phase, disordered b-phase, and metastable c-phase [5–7]. The former two phases have hydrogen bonded zigzag chain structures. The structure of c-phase is in between a- and b-phases. For a-phase, the molecular simulations with potential functions have been performed [2,3,8–10]. These molecular simulation studies have suggested necessity for construction of potential functions based on ab initio potential surfaces.
2. Method The detailed description of the computational procedure of the ab initio MO based lattice energy calculation is given elsewhere [4]. Here, the outline of the method is described. The lattice energy, Elatt , is given as follows: ( ) X X1 eij ; Elatt ¼ ei þ ð1Þ 2 i j6¼i where i runs over all molecules within an Ôasymmetric unitÕ in which the molecule is regarded as an unit object for symmetry operations associated with the space group of crystal. j runs over all the other molecules generated by the symmetry operations. The one-body energy ei is the polarization energy of the molecule i in crystal. The two-body energy eij is the intermolecular interaction energy between molecules i and j. In our method, these energies are directly obtained from ab initio MO calculations.
In this work, the MP2 method is applied to lattice energy calculations. The polarization energy is still evaluated at the HF level and only the pair interaction energy is calculated at the MP2 level using the following equation: n o polðjÞ polðiÞ ei þ ej ; eij ¼ EijMP2 EiMP2 þ EjMP2 ð2Þ where E with the superscript MP2 denotes the total energy of molecule or molecular pair at the MP2 level. Then, the first term of Eq. (2) is the usual intermolecular interaction energy between molecules i and j obtained from supermolecule calcupolðjÞ lations. ei in Eq. (2) is the polarization energy of molecule i involved in the interaction energy between molecules i and j and is calculated at the HF level by the same manner as the one-body polarization energy [4]. 3. Computational details The lattice energy minimization calculations were performed on methanol crystal at the MP2 level with various basis sets; 3-21G, 3-21+G, 3-21G**, 6-31G, 6-31+G, 6-31G*, 6-31+G*, 6-31++G, 6-31G**, and 6-31++G**. We focus on a- and b-phases which have zigzag chains of hydrogen bonds along b-axis in the former and c-axis in the latter (Fig. 1). In the calculations, we assumed ordered structure for both a- and b-phases. Since c-phase is metastable and its structure is unknown [6], no consideration is given to this phase. In order to avoid a multi-minimum problem, the initial parameters including the lattice constants, the molecular orientation, space group, and the number of molecules in unit cell (Z) were taken from the experimental values [5,6] (vide infra). The lattice energy was then minimized with respect to the lattice constants (a, b, and c) and the molecular orientation angles using numerical differentiation. Each methanol molecule in the crystal was kept rigid at the gas phase structure obtained by ab initio MO calculations at the same level. In the lattice energy calculations, we considered from asymmetric unit molemolecules within 20 A cules. The number of molecules within this cut-off length were 664 and 614 for a- and b-phases,
K. Nagayoshi et al. / Chemical Physics Letters 369 (2003) 597–604
599
Fig. 1. The structure of methanol crystals: (a) a-phase and (b) b-phase. The dashed lines show the unit cell.
respectively. Comparing with the lattice energy obtained from the calculation with larger cut-off , the number of methanol molecules: 10 434 for (50 A a-phase and 9606 for b-phase), we estimated that the error due to this approximation is 0.00 and 0.11 kcal/mol for a- and b-phases, respectively, at the MP2/6-31++G** level. Since the error in the lattice energy for a-phase was less than 0.01 kcal/ mol, the Ewald summation [11] was not used. To make even the calculation level and to make sure the difference of phases, the Ewald summation was purposely not used in the calculation of b-phase. The following is the summary of the initial parameters used for the lattice energy minimizations. For a-phase; orthorhombic, P 21 21 21 space group, , and Z ¼ 4 at a ¼ 4:873, b ¼ 4:461, c ¼ 8:867 A 15 K [6]. For b-phase; orthorhombic, Cmcm space , and Z ¼ 4 group, a ¼ 6:43, b ¼ 7:24, c ¼ 4:67 A at 165 K [5]. All calculation were performed with our original program which utilizes the GA U S S I A N 94 [12].
4. Results and discussion For a basis to understand the molecular interaction in the crystal, the optimized structures of
the two isomers of methanol dimer obtained from MP2/6-31++G** calculations are shown in Fig. 2. One is Ôhead to headÕ type (hh) and the other is Ôhead to tailÕ type (ht). The ÔheadÕ means the hydroxyl group and the ÔtailÕ means the methyl group. hh has the linear hydrogen bond between each hydroxyl group and ht has the cyclic hydrogen bond through the two OH groups and the methyl groups. In hh, the O HðOÞ hydrogen and the O–O distance is bond length is 1.888 A . The OH O part in hh is almost linear 2.854 A (\OH O ¼ 172:4°). On the other hand, in ht, there are two weak O HðCÞ type hydrogen bond. The O HðCÞ distance is significantly long ) and the \CH O angle (106.1) is lar(2.401 A gely deviated from 180°. Thus, the Ôhead to tailÕ type of the cyclic hydrogen bonding interaction between two methanol molecules is very weak, even though the hydrogen bond is multiple. For reason of those, the lowest energy configuration of the dimer is linear rather than cyclic. The results are in good agreement with the previous theoretical study [13]. The intermolecular interaction energy of hh ()7.2 kcal/mol) is more stable than that of ht ()3.3 kcal/mol) by 3.9 kcal/mol.This suggests that the Ôhead to headÕ type of hydrogen bonding interaction dominates the crystal packing.
600
K. Nagayoshi et al. / Chemical Physics Letters 369 (2003) 597–604
and Fig. 2. The optimized structures of ðCH3 OHÞ2 at the MP2/6-31++G** level of theory. The bond lengths and angles are in A degrees, respectively. The intermolecular interaction energies, E, are in kcal/mol.
and 11.5°, respectively. The avangles are 0.49 A eraged relative error in the cell lengths is 6.2%. These errors are non-negligibly larger than those of a-phase. This result of b-phase can not be compared with the force field calculations, since they have not been reported. The calculated lattice constants and the molecular orientation highly depend on the basis set employed. We examined the basis set dependence of these quantities by using various basis sets (Fig. 3). The poorest result was obtained at the MP2/3-21G* level; the averaged relative error in the cell lengths exceeded 10% for a-phase. Nevertheless, the order of cell lengths (c > a > b) is well reproduced at most of the basis sets examined. Similar result was also obtained for b-phase.
Table 1 shows the optimized lattice constants of a- and b-phases at the MP2/6-31++G** level. For a-phase, both the lattice constants and the molecular orientation are in good agreement with the experimental values. The maximum deviation of calculated values from the experimental ones [6] is for the lattice constants and 2.2° for the 0.08 A angles of the molecular orientation. The averaged relative error in cell lengths is 0.9%. Note that the results of the force field calculation reported by Mooij et al. [2] are in good agreement with the present results. For b-phase, the optimized lattice constants and the molecular orientation at the MP2/6-31++G** level are in moderate agreement with the experimental values [5]; the maximum deviations of the cell lengths and the orientation
Table 1 Calculated lattice constants and molecular orientation of methanol solid, a- and b-phases, at the MP2/6-31++G** levela a-phase
a b c w h u Lattice energy
b-phase b
c
This work
Force field
Obs
4.85 ()0.02) 4.56 ()0.08) 8.83 ()0.04) 65.4 (+0.2) 85.4 (+2.2) 26.7 (+0.4) )15.0
5.07 (+0.20) 4.66 (+0.02) 8.71 ()0.16) – – – –
4.873 4.641 8.867 65.2 83.2 26.3 –
This work
Obsd
6.08 ()0.35) 7.43 (+0.19) 5.16 (+0.49) 78.5 ()11.5) 185.6 (+5.6) 85.7 ()4.3) )11.4
6.43 7.24 4.67 90.0 180.0 90.0 –
and degrees, respectively. The data in The lattice constants (a; b, and c) and the molecular orientations (w; h, and u) are in A brackets are deduction from observed values. The angles w; h, and u are defined as follows: the rotation by w; h, and u about the molecule lying in the ac plane and C–O axis along the c-axis is rotated to the observed position by w; h, and u about the a-, b-, and c-axes, respectively. b Ref. [2]. The AMMs-II force filed. c Ref. [6]. d Ref. [5]. a
K. Nagayoshi et al. / Chemical Physics Letters 369 (2003) 597–604
Fig. 3. The averaged relative error of the calculated cell lengths. The error is given by fja0 aj=a þ jb0 bj=b þ jc0 cj=cg=3 where a; b, and c are the experimental values and a0 ; b0 , and c0 are calculated ones.
As in the case of a-phase, there is large basis set dependence of the lattice constants, but the order of cell lengths (b > a > c) is again well predicted at most of the basis sets. Although the basis set dependence appears to be not consistent, inclusion of diffuse functions generally reduces the averaged relative error in the cell lengths by 0.5% on average. It is noted that the rather small 321+G basis set gives a good results; the MP2/321+G calculations give the lattice constants in comparable accuracy to the MP2/6-31++G** calculations. Now, let us look at the pair interactions of the molecules in crystal together with their geometries obtained from the MP2/6-31++G** calculation. To simply indicate molecular pairs, they are named as follows; for example, the pair of molecules 1 and 2 (Fig. 1) with the Ôhead to headÕ type interaction in a-phase is called A12-hh, that in b-zphase B12-hh, and so forth. Fig. 4a shows the important molecular pairs in a-phase, whose stabilization energy is more than 1 kcal/mol. Clearly, the most important molecular pair is A12-hh which corresponds to the structure of the most stable dimer (hh) in gas phase. Due to the packing effect, the O–O distance of A12-hh is significantly short
601
) compared to that of the dimer hh in gas (2.697 A ). The other molecular pairs, A13phase (2.854 A ht, A14-ht, and A15-ht, gain rather small stabilization through the distorted Ôhead to tailÕ type interaction. The representative molecular pairs in b-phase are shown in Fig. 4b. A dominant molecular pair is B12-hh which has the Ôhead to headÕ type hydrogen bonding. This pair is resemble to that of A12-hh. As in the case of a-phase, again, one could see the significant shortening of the O–O ) owing to the packing effect, distance (2.740 A is smaller than though the decrease of 0.114 A in a-phase. On the other hand, that of 0.157 A the stabilization energy of the molecular pair B13ht is small due to the weak interaction of Ôhead to tailÕ type. A12-hh and B12-hh form hydrogen bond chain along a certain axis in the crystal. In a-phase, the chain grows along the b-axis (Fig. 5a). The calcu is in reasonable lated O–O distance of 2.697 A . agreement with the experimental value of 2.74 A In b-phase, hydrogen bond chain grows along c-axis (Fig. 5b). The calculated O–O distance of is slightly longer than that in a-phase. On 2.740 A the other hand, the X-ray diffraction study has revealed that the O–O distance of the hydrogen ) is much shorter bond chains in b-phase (2.66 A ). This fact indicates than that of a-phase (2.74 A that there is a significant difference in the packing effect between a- and b-phases. For b-phase, the chain interaction of the Ôhead to headÕ type of hydrogen bonding along c-axis plays a dominant role in the packing. For a-phase, the several weak Ôhead to tailÕ type interactions along a- and c-axes are also important as well as the chain interaction of hydrogen bond along b-axis. Table 2 shows the calculated lattice energy with its component energies and the pair interaction energy of some nearest neighbors. The lattice energy of a-phase ()15.0 kcal/mol) is larger in magnitude than that of b-phase ()11.4 kcal/mol). Although no lattice energy has been observed experimentally, this result is consistent with the fact that a-phase is stable at lower temperature than b-phase. The polarization energy of a-phase is 1.5 kcal/mol larger in magnitude than that of bphase. Similarly, the Coulomb interaction energy
602
K. Nagayoshi et al. / Chemical Physics Letters 369 (2003) 597–604
Fig. 4. The geometrical parameters of molecular pairs in the optimized crystal structure at the MP2/6-31++G** level. The bond and degrees, respectively. The intermolecular interaction energies, E, are in kcal/mol. (a) is a-phase and lengths and angles are in A (b) is b-phase.
of a-phase is 4.4 kcal/mol larger in magnitude than that of b-phase. These results clearly show that the stabilization of a molecule due to surroundings is larger in a-phase than in b-phase. In a-phase, the interaction energy of the dominant molecular pair (A12-hh) provides only 33% of the whole lattice energy. In other word, a large portion of the lattice energy is brought by the many weak interactions. On the other hand, in b-phase, more than 45% of the lattice energy is provided by the dominant
molecular pair (B12-hh). These results also support that there is an apparent difference in nature of the crystal packing between a- and b-phases.
5. Summary We calculated the packing structure of the hydrogen bonding methanol crystal using the ab initio MO based lattice energy minimization
K. Nagayoshi et al. / Chemical Physics Letters 369 (2003) 597–604
603
) in hydrogen bonded chain in the methanol crystal. In the bracket is given the experimental value. Fig. 5. The O–O distance (A (a) a-phase and (b) b-phase.
Table 2 Calculated lattice energy and the intermolecular interaction energies of some nearest neighborsa
Lattice energy Polarization energy Coulomb energy Correlation energy Pair energyb 1-2 1-3 1-4 1-5 1-6 1-7 1-8
a-phase
b-phase
)15.0 )3.5 )14.6 )9.0
)11.4 )2.0 )10.2 )6.7
)4.9 )1.2 )1.1 )1.0 )0.7 )0.5 )0.4
(2) (2) (2) (2) (2) (2) (2)
)5.1 )1.0 )1.0 )0.7 )0.6 )0.2 0.1
(2) (2) (2) (2) (2) (2) (2)
constants, the molecular orientation, and hydrogen bond network in the crystal with satisfactory accuracy. The calculation also revealed the difference in the packing feature between a- and b-phases. Since the examined two phases have different nature of interaction, the simple point charge approximation employed in the estimation of the polarization energy may result in an unbalanced error in the lattice constants. We believe that our computational procedure for the lattice energy would be improved by using the fragment MO method [14,15] which gives the accurate total energy of molecular clusters based on a pair interaction approximation.
a
All energies are given in kcal/mol. Pair energies calculated by Eq. (2). See Fig. 1 for the indices of molecular pairs. The number of the equivalent pairs is given in parenthesis. b
method at the correlated wavefunction level. The calculated lattice constants depended highly on the basis set employed. It has been shown that the MP2/6-31++G** calculations predicts the lattice
Acknowledgements The numerical calculations were partially carried out at the Computer Center of the Institute for Molecular Science. This work was partly supported by ACT-JST.
604
K. Nagayoshi et al. / Chemical Physics Letters 369 (2003) 597–604
References [1] John Maddox, Nature 335 (1988) 201. [2] W.T.M. Mooij, B.P. van Eijck, J. Kroon, J. Phys. Chem. A. 103 (1999) 9883. [3] M. Hloucha, A.K. Sum, S.I. Sandler, J. Chem. Phys. 113 (2000) 5401. [4] T. Ikeda, K. Nagayoshi, K. Kitaura, Chem. Phys. Lett., in press. [5] K.J. Tauer;, W.N. Lipscomb, Acta Cryst. 5 (1952) 606. [6] B.H. Torrie, S.-X. Weng, B.M. powell, Mol. Phys. 67 (1989) 575. [7] B.H. Torrie, O.S. binbrek, M. Strauss, I.P. Swainson, J. Solid State Chem. 166 (2002) 415. [8] A.K. Sum, S.I. Sandler, R. Bukowski, K. Szalewicz, J. Chem. Phys. 116 (2002) 7627.
[9] R. Bukowski, K. Szalewicz, C.F. Chabalowski, J. Phys. Chem. A. 103 (1999) 7322. [10] J. Alonso, F.J. Bermejo, M. Garcia-Hernandez, J.L. Martines, W.S. Howells, A. Criado, J. Chem. Phys. 96 (1992) 7696. [11] P.P. Ewald, Ann. Physik. 64 (1921) 253. [12] M.J. Frisch et al., GA U S S I A N 94, Revision D.4, Gaussian, Pittsburgh, PA, 1995. [13] R.A. Provencal, J.B. Paul, K. Roth, C. Chapo, R.N. Casaes, R.J. Saykally, G.S. Tschumper, H.F. Schaefer III, J. Chem. Phys. 110 (1999) 4258. [14] K. Kitaura, T. Sawai, T. Asada, T. Nakano, M. Uebayashi, Chem. Phys. Lett. 312 (1999) 319. [15] K. Kitaura, E. Ikeo, T. Asada, T. Nakano, M. Uebayash, Chem. Phys. Lett. 313 (1999) 701.