Liquid methanol Monte Carlo simulations with a refined potential which includes polarizability, nonadditivity, and intramolecular relaxation

Liquid methanol Monte Carlo simulations with a refined potential which includes polarizability, nonadditivity, and intramolecular relaxation Maximiliano Valde´z-Gonza´lez,aHumberto Saint-Martin, and Jorge Herna´ndez-Cobos Instituto de Ciencias F´sicas, Universidad Nacional Auto´noma de Me´xico, Apartado Postal 48-3, 62251 Cuernavaca, Morelos, Mexico Regla Ayala and Enrique Sanchez-Marcos Departamento de Qu´mica F´sica, Universidad de Sevilla, 41012-Sevilla, Spain Ivan Ortega-Blakeb Departamento de F´sica Aplicada, Cinvestav, Km. 6 Antigua Carretera a Progreso, Cordemex, Me´rida 97310, Yucata´n, Mexico Received 9 March 2007; accepted 2 October 2007; published online 12 December 2007 Monte Carlo simulations of liquid methanol were performed using a refined ab initio derived potential which includes polarizability, nonadditivity, and intramolecular relaxation. The results present good agreement between the energetic and structural properties predicted by the model and those predicted by ab initio calculations of methanol clusters and experimental values of gas and condensed phases. The molecular level picture of methanol shows the existence of both rings and linear polymers in the methanol liquid phase. © 2007 American Institute of Physics. DOI: 10.1063/1.2801538 I. INTRODUCTION Nowadays it is possible to use refined potentials in numerical simulations of physicochemical systems involving small molecules1–3and attain sufficiently good agreement with the experimental observations. The reliability that such validation confers upon the simulations allows for the construction of a molecular image that contributes to a better understanding of the phenomena. It seems that construction of refined potentials requires paying attention to molecular properties such as polarizability, intramolecular relaxation, and nonadditivity.4–7These flexible potentials can be adjusted to ab initio data of the molecule and the intermolecular interaction, with no reference to a particular thermodynamic state of the condensed phase, allowing for its unbiased use in any condition. A great effort has been made for the development of water potentials1,4,8–10 and in some cases for other small molecules.11,12 It is indeed convenient to extend the use of this tool to other systems and also to determine to what extent the inclusion of different molecular properties that add dearly to the computational cost is required for the proper reproduction of the experimental data.13 Liquid methanol is of great interest given its many uses, particularly that as a common organic solvent and more recently as an important fuel alternative.14 Additionally, a lot of work has been devoted to the understanding of its unusual physical properties15 that have been associated with a peculiar molecular behavior that is conducive to methanol being certainly one of the most structured liquids. In the crystal phase methanol presents long one-dimensional chains of hydrogen bonds.16,17 It has therefore been thought that something similar could be occurring in the liquid phase and be responsible for its peculiar behavior. Several numerical simulations18–23 and neutron diffraction experiments24–26 have produced data that support this view. On the other hand, experimental results of neutron diffraction,27,28 x-ray scattering,29,30 and x-ray emission spectroscopy31 have been taken to support the existence of cyclic clusters of methanol in the liquid phase. Kashtanov et al.31 have suggested that numerical simulations could be using potentials that are not able to reproduce the hydrogen bonding network in ring structures. There is a substantial number of numerical simulations of methanol18–23,32–34 and some of them which include polarizability35–38 have shown the need for refined potentials. In this work a refined methanol-methanol potential that uses the mobile charge densities in harmonic oscillators MCDHO model4which includes polarizability, nonadditivity, and intramolecular relaxation is presented. The potential is adjusted to ab initio surfaces and tested in the reproduction of ab initio methanol clusters, their energies, and structures. The potential is then used in Monte Carlo simulations of liquid methanol. This serves to test the potential and determine the validity of its use. In addition we tested how accounting for different molecular properties, such as polarizability and intramolecular relaxation, affects the reproduction of the experimental properties. The quality of the potential helps us to further the understanding of the behavior of liquid methanol. aElectronic mail: [email protected] bOn leave from Instituto de Ciencias Fı´sicas, Universidad Nacional Auto´nomadeMe´xico. THE JOURNAL OF CHEMICAL PHYSICS 127, 224507 2007 0021-9606/2007/127 22 /224507/14/$23.00 © 2007 American Institute of Physics127, 224507-1 R eu s e of AIP Pub lis h i ng c ontent is s ub j e c t to the te rms : http s ://pub lis h i ng.a i p.o r g/autho rs / ri ght sand - pe rmissi on s . D o w n l oaded to IP: 150.214.182.116 On: W ed, 19 O c t 2016 14:19:58 II. METHODS A. Potential energy surfaces The methanol-methanol potential was developed by reproducing ab initio energy surfaces for the intramolecular relaxation, the dipole moments of different methanol structures, the pairwise interaction, and the nonadditive terms of the intermolecular interaction. Molecular orbital calculations for all cases were done with the GAUSSIAN98 program.39 We considered in the calibration procedure Gaussian basis sets 6-31+ +G** and 6-311+ +G** and correlation-consistent basis sets cc-pVDZ, aug-cc-pVDZ, cc-pVTZ, and cc-pVQZ.40–42 Of course a compromise has to be made between the computational cost and the quality of the molecular calculation due to the fact that a great number of structures need to be considered. Hence, in order to determine the optimal level of calculation we looked into the prediction of various properties reported at different levels of molecular theory. Table I shows the optimal structures for the monomer of methanol predicted with different basis set sizes at the MP2 and QCC levels. It is clear that there is a fast convergence of the structural parameters and that the aug-cc-pVDZ basis set yields a good approximation to more expensive and complete basis sets. In a previous study Wang et al. compared the vibrational spectra of methanol predicted with different basis sets.43 It was found that a basis set of similar size yielded a reasonable approximation to the experimental spectra. The dipole moment shows more discrepancy, with a value 10% larger for the most extended basis set. However, since the value predicted by the aug-cc-pVDZ basis set is closer to the experimental dipole moment44 of 1.69 D, we decided to use this basis set for the description of the intramolecular surface. Table II shows the structural parameters of the optimal methanol dimer predicted by different basis sets. It is clear that there is convergence in the predicted dimer configuration. The corresponding energies are also presented and compared to other data in the literature. The values corresponding to the three largest basis sets were computed with structures optimized with a 6-31+ +G** basis set and monomers kept frozen in the dimer optimization.45 Even if they are not fully comparable they provide a good reference point on the expected convergence value, and since they correspond to a partial optimization it is clearly a lower limit. The optimal dimers for the six smaller basis sets were optimized in their own basis set and counterpoise CP correction apTABLE I. Structural parameters for the minimum energy methanol monomer predicted at the MP2 and QCC levels. Distances are in angstroms, bond and dihedral angles dare in degrees, and dipole moment in Debye. Method r CO rOH COH dHCOH MP2/cc-pvDZ 1.4171 0.9656 106.3 179.88 1.63 MP2/6-31+ +G** 1.4290 0.9644 108.5 179.88 1.99 MP2/6-311+ +G** 1.4217 0.9544 107.3 179.88 1.93 MP2/aug-cc-pVDZ 1.4342 0.9659 107.9 180.00 1.71 MP2/cc-pvTZ 1.4188 0.9594 107.4 179.83 1.65 MP2/aug-cc-pvTZ 1.4239 0.9611 108.0 179.81 1.70 MP2/cc-pVQZ 1.4177 0.9577 108.0 179.83 1.68 MP2/aug-cc-pVQZ 1.4353 0.9658 107.9 179.6 1.70 QCC TZV 2p,2d ++a1.4286 0.9583 107.6 180.0 1.81 aReference 34. TABLE II. Structural parameters of the optimal methanol dimer predicted at the MP2 level left columns . Distances rare in angstroms, and bond and dihedral angles dare in degrees. Predicted interaction energy kcal/mol for the optimal methanol dimer at the MP2 and CCSD Tlimit levels right columns .E int corresponds to the energy of the dimer minus the energy of the relaxed monomers, Edef is the deformation energy, and Eint =Eint+Edef is the interaction energy. The relative time for a single point calculation for each basis set is also presented. Basis set r O··H r O··O OH··O HO··H d COH··O d OH··OC Eint Edef Eint Rel. time cc-pVDZ 1.887 2.819 159.3 101.7 109.1 −20.7 −3.53 0.20 −3.33 1.0 6-31+ +G** 1.887 2.853 172.5 117.9 109.0 16.9 −5.22 0.11 −5.11 1.5 6-311+ +G** 1.886 2.846 171.9 122.5 98.3 30.0 −4.97 0.11 −4.86 5.9 aug-cc-pVDZ 1.887 2.847 168.2 112.7 132.3 −6.7 −5.22 0.11 −5.11 15.7 cc-pVTZ 1.872 2.802 160.4 131.0 97.8 24.9 −5.09 0.19 −4.90 39.0 aug-cc-pVTZ 1.877 2.836 169.4 117.0 129.9 1.0 −5.60 0.11 −5.49 3486.3 cc-pVQZ −5.21a¯ cc-pV5Z −5.39a¯ CCSD Tlimit b−5.45a¯ aReference 45. baug-cc-pVTZ optimized geometry. 224507-2 Valde´z-Gonza´lez et al. J. Chem. Phys. 127, 224507 2007 R eu s e of AIP Pub lis h i ng c ontent is s ub j e c t to the te rms : http s ://pub lis h i ng.a i p.o r g/autho rs / ri ght sand - pe rmissi on s . D o w n l oaded to IP: 150.214.182.116 On: W ed, 19 O c t 2016 14:19:58 plied to the energy obtained by subtracting the relaxed methanol energies. The convergence limit has also been validated by comparison to the experimentally derived binding energy.46 With this in mind, we can estimate that the value obtained by the aug-cc-pVDZ basis set has a small underestimation of 5% and considering the relative computational costs presented in Table II, the aug-cc-pVDZ basis set was chosen for the determination of the dimer interaction energy at the MP2 CP corrected level. Since CP correction applied to fully optimized dimers including monomer relaxation is not well defined, we performed the following algorithm for the estimation of CP correction at the aug-cc-pVDZ basis set level. The dimer interaction was computed by subtracting the energies of the resulting methanol monomers in the fully optimized dimer from the energy of the dimer. The energies of these monomers were computed in the complete basis of the dimer. Then the deformation energy of each of these monomers with respect to the optimal one was computed in the same basis set and subtracted from the previous energy; using a much larger basis set as the aug-cc-pVQZ did not produce any significant difference in the deformation energy values. This algorithm allows for a better search of the optimal structure that does not necessarily correspond to the minimum of the energy surface computed ab initio, due to CP correction. Of course it is possible that CP overcorrects the interaction energy, but since CP is in itself small, this effect is negligible. The monomer deformation energy i, the dimer interactionenergyV ij, and the three-body nonadditivity ijk was calculated using the many-body expansion described in Ref. 47. The three-body nonadditivities were computed with the SCF CP /aug-cc-pVDZ level. The reason why correlation energy was not considered in these calculations is that it has been shown that correlation energy is quite additive3,48 and restriction to the uncorrelated level entails significantly less computer time. The same level of calculation was used to compute some four-body nonadditive terms in order to assess the magnitude of these contributions. They were found to be small enough—and reasonably well reproduced by the potential—to not merit the effort of adjusting to the whole four-body potential energy surface. The sample points in the potential energy surfaces were chosen in an iterative manner. Thus, initial samples were taken from a regular scan of the hypersurface with respect to the different degrees of freedom. An initial potential parametrization was fitted to the preliminary samples and employed to predict, via numerical simulations at standard temperature, more monomers, dimers, and trimers in the gas phase as well as clusters appearing in the liquid phase. These structures were then computed at the ab initio level and added to the potential energy surface. This procedure continued until self-consistency was attained with the same accuracy as that obtained in the previous fitting. In this way, for instance, the optimal monomer energy predicted by the potential differs in only 0.06 kcal/mol with respect to that computed at the ab initio level. The final surfaces to be fitted consisted of athe polarizability and dipole moment of the optimal monomer, b555 monomer structures where both the deformation energy, relative to the optimal monomer, and the dipole were considered for fitting, c762 dimer structures where the interaction energy was fitted, and d153 trimer structures where the three-body nonadditive contributions where fitted. The vibrational spectra presented in this work were calculated using the normal mode theory.49 In the condensed phase, the semiclassical method suggested by Reimers and Watts for water and ice50 was implemented. The spectra was calculated by averaging over 200 configurations taken from a Monte Carlo calculation with 500 molecules in the simulation cell. The configurations are separated by 100 000 Monte Carlo steps. B. Model potential The model considered in this work uses a functional form based on the MCDHO model,4as it allows for the inclusion of intramolecular flexibility, nonadditivity, and polarizability. The electron cloud of each atom of a molecule is represented as a negative mobile charge density with radial exponential decay, attached to a positively charged point by a harmonic oscillator to simulate the interaction between charges of atoms forming a chemical bond. All other chargecharge interactions consider all charges as points instead of densities. The specific expression for the water-water interaction was given in Eqs. 4 – 6 of Ref. 4.Herewegiveamore general formulation, suitable to be used with larger molecules. In the following equations rij= r i−rjis the distance between the corresponding centers, either fixed positive nuclei Z or mobile negative charges q, where i and j subscripts correspond to sites i and j. The intramolecular energy is composed of the following terms. 1The electrostatic interaction between atomic centers: UZi,Zj=ZiZj rij . 1 2The electrostatic interaction between mobile charge densities qiand the atomic center charge Zjof the bonded atoms, where iis the intramolecular decay length of the charge density: Uqi,Zj=qiZj rij 1−rj i +1 exp −2rij i . 2 The mobile charge densities are considered as point charges when interacting with any other atomic center. In these cases the electrostatic interaction is given by Uqi,Zj=qiZj rij . 3 3The interaction between mobile charges with decay lengths iand jattached to bonded atoms. In this case a two-center integral should be used; however, an approximate expression that gives good accuracy at all relevant distances was found, i.e., 224507-3 Liquid methanol Monte Carlo simulations J. Chem. Phys. 127, 224507 2007 R eu s e of AIP Pub lis h i ng c ontent is s ub j e c t to the te rms : http s ://pub lis h i ng.a i p.o r g/autho rs / ri ght sand - pe rmissi on s . D o w n l oaded to IP: 150.214.182.116 On: W ed, 19 O c t 2016 14:19:58 U qi,qj=qiqj rij 1−i j +i+j2 i+j3rij +1 exp −2i j +i+j2 i+j3rij , 4 again the mobile charge densities are considered as point charges when interacting with mobile charge densities of nonbonded atoms Uqi,qj=qiqj rij . 5 4A Morse potential between each pair of bonded atoms, with depth Dij, inverse decay length ij, and equilibrium parameter rij eq, Ubri,rj=D ij exp −2ij rij −rij eq −exp −ij rij −rij eq . 6 5A quadratic plus Urey-Bradley UB terms for each bond angle ijk, with parameters k and kUB, and the equilibrium values ijk 0and rik 0, U=k 2ijk −ijk 02, 7 UUB =kUB 2rik −rik 02. 8 6A periodic function over the dihedral angle between four consecutive bonded centers Udih =k 21+cos n + . 9 7Exponential terms between the nonbonded atoms, i.e., that do not form angles or covalent bonds. Unb =Aexp −ar +Bexp −br .10 Hence the analytical expression to compute the intramolecular energy is US= i S j i U Zi,Zj+U qi,qj+ i S j i U qi,Zj + i S 1 2kirii 2+ bonds UbZi,Zj+ angles U+U UB + dih Udih + nonbond ULJ,11 where the third term corresponds to the harmonic energy of the mobile charge at a distance rii from its nucleus and the summation i S refers to the entire molecule. TABLE III. Parameters of the model potential. Distances are in a.u., angles in deg, and energy is in kcal/mol. The different parameters are described in the text. Rcstands for the hard core cutoff radius used during simulation and are in a.u. Site Site Site Site Z q k C 3.900 283 −3.408 234 569.650 1.223 331 O 0.366 121 −1.398 671 435.404 1.597 584 HC0.596 934 −0.601 763 602.237 0.900 642 HO1.990 457 −1.442 636 455.769 0.849 260 Dij ij rij eq O C 4.595 276 0.977 601 3.997 212 HCC 70.427 959 1.078 138 2.292 126 OH O21.473 454 1.584 645 2.021 717 kijk 0kUB rik 0 HCC O 89.700 109.471 7.957 3.955 COH O42.007 107.800 26.103 3.782 HCCH C55.768 109.000 14.745 4.063 kn HCCOH O0.118 728 3 0.0 Aij aij Bij bij Rc HCHO480.231 97 0.176 31 −892.495 99 0.047 74 C C 97 519.618 11 2.094 63 −2.294 13 0.266 26 5.6 C O 6 496.790 10 1.276 07 −1782.699 82 0.956 43 5.3 CH C1 644.359 38 1.721 79 10.903 08 0.394 32 4.2 CH O435.951 37 0.985 75 −2.765 14 0.119 54 3.7 O O 40 041.925 58 1.065 55 −36 415.261 63 1.047 97 4.3 OH C1 634.413 65 1.477 24 0.137 03 0.024 35 3.7 OH O4 001.505 20 2.034 11 8.159 27 0.299 31 2.8 HCHC2 904.622 84 2.096 89 −9.006 16 0.478 64 2.83 HCHO1 114.574 06 1.201 24 −1132.260 23 1.178 35 2.83 HOHO623.490 48 1.193 27 −347.755 15 0.965 90 2.95 224507-4 Valde´z-Gonza´lez et al. J. Chem. Phys. 127, 224507 2007 R eu s e of AIP Pub lis h i ng c ontent is s ub j e c t to the te rms : http s ://pub lis h i ng.a i p.o r g/autho rs / ri ght sand - pe rmissi on s . D o w n l oaded to IP: 150.214.182.116 On: W ed, 19 O c t 2016 14:19:58 The intermolecular energy is composed of the following terms. 1The electrostatic interactions between atomic centers plus exponential terms with parameters Aij,B ij,a ij,and bij; Uinter Zi,Zj=A ij exp −aijrij +B ij exp −bijrij +ZiZj rij . 12 2The electrostatic interaction between mobile charges qi considered as point charges and charges Zj, Uinter qi,Zj=qiZj rij .13 3The electrostatic interaction between mobile charges, Uinter qi,qj=qiqj rij .14 Therefore, the energy of a cluster with N molecules is given by U= S=1 N T=1 S−1 i S j T Uinter Zi,Zj+U inter qi,Zj +U inter qj,Zi+U inter qi,qj+U S,15 where summations i S and j T refer to entire molecules. The interaction energy requires the subtraction of the intramolecular energies of the isolated molecules, US 0, U=U− S=1 N US 0,16 thus taking into consideration the energetic cost of polarizing and deforming each molecule in the cluster or the condensed phase. Unlike the case of water,4and because of the complexity of the potential energy surface due to the increased number of degrees of freedom, it was found convenient to use hard core cutoff radii for the interactions as part of the potential. So in addition to the set of parameters fitted to reproduce the surfaces, using the program VA 05AD ,51 the potential definition TABLE IV. Structural parameters for the optimal methanol monomer. Comparison of MCDHO results against ab initio and experimental results. Error bars for the MCDHO results at 298 K correspond to 2 using the method of Flyvbjerg and Petersen, Ref. 71, for the computation of . Distances are in angstroms, bond and dihedral angles dare in deg, and dipole moment in Debye. Method r CO rOH COH dHCOH aug-cc-pVDZ 1.4342 0.9659 107.9 180.00 1.71 MCDHO 0K 1.4386 0.9639 107.8 180.00 1.72 MCDHO 298 K 1.4423±0.0014 0.9667±0.0004 108.4±0.3 ¯1.72±0.01 Expt.a1.434 0.937 105.93 ¯1.69 Expt.b1.4246±0.0024 0.9451±0.0034 108.53±0.48 ¯¯ aReference 44. bReference 72. FIG. 1. Comparison between the monomer deformation energies , monomer dipole moments ,dimer interaction energies as described in Ref. 47 Vij , and three-body nonadditivities as described also in Ref. 47 ijk predicted at the ab initio level and those predicted by the model for the same structures. 224507-5 Liquid methanol Monte Carlo simulations J. Chem. Phys. 127, 224507 2007 R eu s e of AIP Pub lis h i ng c ontent is s ub j e c t to the te rms : http s ://pub lis h i ng.a i p.o r g/autho rs / ri ght sand - pe rmissi on s . D o w n l oaded to IP: 150.214.182.116 On: W ed, 19 O c t 2016 14:19:58 includes a set of cutoff values. Differentiating the methyl hydrogen from the hydroxyl hydrogen also proved convenient, so they were treated as different atomic species. III. RESULTS A. Model potential performance 1. Reproduction of the adjusted properties The fitted parameters for the methanol-methanol potential and the corresponding cutoff values are presented in Table III. The quality of the reproduction of the monomer deformation energy surface and dipole moment is presented in Fig. 1, as well as the two-body interactions and the threebody nonadditivity. In addition the trace of the polarizability tensor also fitted by the model comes to be xx=21.2, yy =23.3, and zz=18.8 a.u. which compares rather well to the ab initio values xx =20.5, yy=23.2, and zz =19.8 a.u. and the experimental isotropic values of of 21.8 Ref. 52 and 22.0 a.u.53 2. Reproduction of molecular clusters A comparison between the geometries of the optimal methanol monomer predicted by the model, the corresponding ab initio monomer, and the experimental values is presented in Table IV. A very good agreement of the model with the optimal ab initio structure was found. Furthermore, the Monte Carlo simulation of the single monomer at 298 K predicts geometrical values which are very close to the experimental ones. In Fig. 2the vibrational spectra of the methanol molecule in the liquid phase is presented and compared to the experimental data. It was computed via a normal mode analysis. In order to estimate the reliability of this treatment we compared its predictions for the water molecule with those estimated from the velocity autocorrelation in a molecular dynamics MD simulation with the MCDHO potential by Stern et al.54 We can see in Table Vthat the two methods yield similar results for the shift from the gas to liquid phase, giving confidence in the application of the normal mode analysis to the methanol molecule. In Table VI the vibrational spectra predicted for the methanol molecule in TABLE V. Vibrational frequencies of water in gas and liquid phases. All frequencies are in cm−1. Expt. MCDHO MD a MCDHO NM b Expt. Liquid-gas MCDHO MD Liquid-gas MCDHO NM Liquid-gasGascLiquiddGas Liquid Gas Liquid OH symm. stretch 3756 3918 3668 3928 3742 −250 −186 3557 −150 OH antisymm. stretch 3657 3785 3407 3789 3346 −378 −443 HOH bend 1595 1670 1651 1740 1660 1754 75 89 94 aReference 54. bNormal modes. This work. cReference 73. dReference 74. FIG. 2. Vibrational spectral density of liquid methanol, MCDHO, and experimental results. 224507-6 Valde´z-Gonza´lez et al. J. Chem. Phys. 127, 224507 2007 R eu s e of AIP Pub lis h i ng c ontent is s ub j e c t to the te rms : http s ://pub lis h i ng.a i p.o r g/autho rs / ri ght sand - pe rmissi on s . D o w n l oaded to IP: 150.214.182.116 On: W ed, 19 O c t 2016 14:19:58 the gas and liquid phases are compared against the ab initio and experimental values. It can be seen that there is a reasonable agreement in the gas phase description, but less so for the liquid phase. In particular, the shift from gas to liquid phase for the OH bond stretch mode is much reduced. The optimal methanol dimer predicted by the model is in very good agreement with the optimal ab initio dimer, as shown in Table VII.InFig.3inset a superposition of ab initio and model optimum dimers is presented showing a discrepancy in a rotation along the hydrogen bond axis, something very difficult to prevent, e.g., it also occurred in the water dimer.4A search for an ab initio CP corrected dimer as described previously produced the parameters presented in the third column of Table VII, which show a better agreement to the experimental values. There is also an improved agreement between the model optimal dimer and the CP corrected dimer. Of course it is important to reproduce the hydrogen bond interaction not only in the single optimal structure. In Fig. 3 the interaction energy profile for the optimal approach of the methanol monomers along the hydrogen bond is presented showing an excellent agreement between the model and the ab initio curve. In Table VIII and Fig. 4we compare the optimal methanol trimers predicted by the model and the ab initio calculations. In this case it was not possible to find the CP corrected optimal structure since the ab initio optimization follows the CP uncorrected surface. But the optimal trimer predicted by the model when calculated at the ab initio level with CP correction turned out to have a mere 0.14 kcal/mol difference with the one obtained by the ab initio minimization. Hence trimer structure reproduction is in good agreement with the ab initio calculation except for a small overestimation of the interaction energy and an overestimation of 0.05 A ˚in the oxygen-oxygen distance. It is also important to look into the hexamer ring, since there is a discussion regarding the importance of this structure in the liquid phase of methanol, as well as a report31 on the possibility that simple potentials do not reproduce well this structure. In Table VIII we compare some structural parameters predicted by the model and ab initio calculations at the MP2/6-31+ +G** level for the optimal hexamer. The table shows a general agreement between both structures. There is an apparent elongation of the hydrogen bond distance in the model predicted structure, but this elongation of 0.04 A ˚is quite similar to the apparent elongation between the CP corrected and uncorrected structures at the ab initio level. Furthermore, the level of ab initio calculation is smaller than the one used for computing the energy surface, because of the computational cost involved, hence the agreement is quite adequate. From this it can be concluded that the model potential is reproducing the methanol clusters accurately and that it can be used to study the condensed phase. B. Numerical simulations The results obtained in reproducing the molecular properties, the interaction energies, and structures of small clusters support the reliability of the model potential in studies of the condensed phase. Liquid methanol simulations were carried out with the Monte Carlo Metropolis algorithm in the canonical ensemble NVT with periodic boundary conditions and Ewald summation. A cubic box of length 2.3437 A ˚ with 500 molecules was used to reproduce a standard density of 0.782 g/cm3at 298.15 K. The procedure of updating only the polarization of the trial molecule called single update has been criticized for lacking the condition of detailed balance55,56 that is sufficient but not necessary for a valid sampling;57 however, apart from hindering the convergence of the dipole-dipole correlation function,58 from which the dielectric constant can be computed, the single update scheme does not produce any significant error compared to algorithms that comply with detailed balance,55,56,59,60 but substantially increase the computational cost. Thus single update has been used in this work. 100 106configurations were used to attain equilibrium starting from a randomly generated configuration, and 100 106configurations were used for extracting the average values. We checked that both intramolecular and intermolecular degrees of freedom were equilibrated. For this we considered separately each energy for 50 106configurations and both energies were at equilibrium. TABLE VI. Vibrational frequencies of methanol in gas and liquid phases. All frequencies are in cm−1. aug-cc-pVDZ Expt.aMCDHO Expt. Liquid-gas MCDHO Liquid-gasGas Gas Liquid Gas Liquid O–H stretch 3839 3681 3328 3752 3679 −353 −73 d-stretch 3189 3000 2980 3540 3478 −20 −62 3130 2960 2946 3260 3276 −14 16 CH3s-stretch 3052 2844 2834 3252 3101 −10 −151 CH3d-deform 1505 1477 1480 1592 1617 3 25 1494 1477 1480 1538 1617 3 79 CH3s-deform 1466 1445 1450 1521 1557 5 36 O–H bend 1367 1345 1418 1240 1492 73 252 CH3d-rock 1170 1165 1163 1167 1185 −218 1076 1065 1115 1116 1093 50 −23 C–O stretch 1046 1033 1030 981 973 −3−8 Torsion 315 200–295 655 307 658 455–360 351 aReference 75. 224507-7 Liquid methanol Monte Carlo simulations J. Chem. Phys. 127, 224507 2007 R eu s e of AIP Pub lis h i ng c ontent is s ub j e c t to the te rms : http s ://pub lis h i ng.a i p.o r g/autho rs / ri ght sand - pe rmissi on s . D o w n l oaded to IP: 150.214.182.116 On: W ed, 19 O c t 2016 14:19:58 The enthalpy of vaporization was computed as Hvap =U vap −Uliq +RT, 17 where Uvap is the energy of the vapor phase and was calculated from a Monte Carlo simulation with a single molecule at 298 K in the same cubic box used for the liquid and using 1106configurations for equilibration and 5 106configurations to average the energy of the system. The vaporization enthalpy predicted Hvap=8.81±0.04 kcal/mol, which compares rather well with the reported experimental value of Hvap=8.946±0.005 kcal/mol.61–63 The predicted molecular properties of the methanol molecule in the liquid phase and their experimental counterparts are presented in Table IX. There is a general agreement between the predicted and the experimental values except for the dipole moment which appears to be underestimated. However, there is an ample discussion in the literature suggesting that the experimental value of 2.85 D is overestimated. Pieruccini and Saija64 considered that the dielectric constant of 33 indicates a dipole in the liquid of 2.39 D. The same value has been proposed by Wick and Dang65 in a classical simulation and Handgraff and Meijer66 in a Car and Parinello simulation. Similarly Weerasinghe and Smith.67 with an empirical model and Martı´n et al.38 in a quantum mechanics/molecular mechanic simulation propose a value of 2.4 D. In summary the value predicted by the model is quite good and therefore all properties of the molecule in the liquid are well reproduced. The total radial distribution function predicted by the model was estimated and compared with those of Yamaguchi et al.25 and Adya et al.26 in Fig. 5.Thefirst two curves were constructed from the partial radial distribution functions using the weighting scheme suggested by Adya et al.26 In Fig. 5the comparison shows that the disagreement between the TABLE VII. Structural parameters for the minimum energy methanol dimer. Comparison of MCDHO results against ab initio results. The third column corresponds to the counterpoise corrected dimer found as described in text. Distances are in A ˚, bond and dihedral angles dare in deg, and energy in kcal/mol. MCDHO MP2/aug-cc-pVDZ MP2/aug-cc-pVDZ Counterpoise corrected Expt. rO··H 1.940 1.887 1.967 1.96±0.02a rO··O 2.892 2.847 2.920 ¯ OH¯O166.0 168.2 166.2 ¯ HO¯H101.5 112.7 112.4 ¯ dCOH··O 159.8 132.3 131.3 ¯ dOH··OC −53.6 −6.7 −6.4 ¯ U−5.27 −5.11 −5.15 ¯ aReference 76. FIG. 3. Methanol-methanol dimer interaction energy along the hydrogen bond approaching line. MCDHO predicted values stars vs ab initio results circles . Inset: optimal methanol dimer predicted by the model top right and the ab initio top left one. The geometries were superimposed on the oxygen atom and as close as possible to the oxydryl hydrogen and the carbon atom of the bottom monomer. 224507-8 Valde´z-Gonza´lez et al. J. Chem. Phys. 127, 224507 2007 R eu s e of AIP Pub lis h i ng c ontent is s ub j e c t to the te rms : http s ://pub lis h i ng.a i p.o r g/autho rs / ri ght sand - pe rmissi on s . D o w n l oaded to IP: 150.214.182.116 On: W ed, 19 O c t 2016 14:19:58 theoretical and the experimental lines is similar to that shown between both experimental curves. This supports the validity of the potential; even if the theoretical peak is still slightly shifted towards longer distances. We can go further and look into the radial distribution functions for different atomic pairs and compare to the experimentally derived curves of Yamaguchi et al.25 Figs. 6 and 7. There is a general agreement of the model radial distribution functions to those obtained from the experiment, except for the oxygen-oxygen radial distribution function, the one involved in the hydrogen bonding between methanol molecules, where the model presents a firstpeakatadistance about 0.1 A ˚longer compared with the curves presented by Yamaguchi et al.,25 and 0.08 A ˚longer compared to the curve derived from the data presented by Adya et al.26 There is also a marked difference in the first minimum and the second peak, the experimental curve being more structured. These discrepancies are rather surprising considering the good agreement with all previous comparisons, either with ab initio or experimental counterparts. The possibility that the potential was not properly reproducing the nonadditivity in long hydrogen bonded structures, as suggested by Kashtanov et al.,31 was considered, but it was shown that hexamers were properly reproduced. We also considered the fact that nonadditivity at the MP2 level was not computed because it is known that correlation energy is highly additive,48 but the possibility remains that, in these cases, these corrections were not negligible. Xantheas68 reported a substantial effect of MP2 corrections on water nonadditivity, but in that work a comparison is made between the optimal structures of the trimer predicted at the MP2 and those predicted at the selfconsistent field SCF level; indeed, pair interaction and thus MP2 is crucial for structure determination, leading to different structures with different nanoadditivity values. Nonetheless in order to check if this could be a source of error, we computed the three-body nonadditivity surface at the MP2 level and compared it with the corresponding SCF surface. No meaningful difference appeared at any point. We also considered the possibility that the cell size was not large enough and it was hindering the formation of larger clusters. However, the radial distribution functions obtained with a cell containing 1000 molecules rendered the same results as those from the cell containing 500 molecules. It is clear that discrepancies between theoretical and experimental radial distribution functions are related to the hydrogen bonding. In a recent work it has been shown that both the radial distribution functions and the dipole moment are quite sensitive to pressure.65 Hence the observed differences could be due to a diminished pressure in the NVT simulation at the fixed experimental density. In order to check this a NPT simulation of the same system was performed at 1 atm. The results show that the density obtained increases 0.03 g/cm3with respect to the previous experimental value of 0.782 g/cm3, and the vaporization enthalpy becomes Hvap =8.99 kcal/mol, getting closer to the experimental value. But the radial distribution functions did not change in any appreciable manner. The only remaining explanation would be that since the level of ab initio calculations is slightly underestimating the binding energy, this would conduce to less compact structures. We have to recall that a small shift for the OH bond stretch mode was predicted between the gas and liquid phases, a deficiency that could be involved in the discrepancy observed. However, it is surprising that the vaporization energy is well reproduced being a very sensitive parameter. It could also be that the experimental biatomic radial distribution functions are reflecting some model dependency since the empirical potential structure refinement approach was used.25 The model oxygen-oxygen radial distribution function indicates well defined first neighbors, but less defined second neighbors when compared to the equivalent water radial disTABLE VIII. Structural parameters for the optimal methanol trimer. Comparison between MCDHO and ab initio results at the MP2/aug-cc-pVDZ counterpoise corrected level left columns and geometrical parameters for the methanol hexamers right columns . Ab initio MP2/6-31+ +G** and MCDHO results. Oaand Odstand for the acceptor and donor oxygen in the hydrogen bonds. Distances are in A ˚, bond and dihedral angles d are in deg, and energies in kcal/mol. The model values correspond to the average and standard deviation values produced by slightly different geometries observed in the hexamer. Trimer Ring hexamer Chain hexamer MCDHO Ab initio MCDHO Ab initio MCDHO Ab initio rO··H 11.909 1.853 r O–O 2.72± 0.01 2.68 rO··H 21.959 1.881 O–O–O 109.5±4.1 117.5 rO··H 31.889 1.860 Oa–Od–C 108.2± 1.7 106.3 rO··O 12.819 2.766 Od–Oa–C 113.5± 3.6 119.4 rO··O 22.853 2.790 d O–O–O–O 59.7±6.5 30.7 rO··O 32.800 2.766 d C–O–O–C 167.7±6.2 123.1 C–O··H 1120.90 127.303 C–O··H 2113.14 109.864 C–O··H 3113.29 109.864 dO–H··O–H 1−6.51 −14.022 dO–H··O–H 2−4.48 −7.972 dO–H··O–H 32.81 −7.888 U−16.18 −15.22 −42.89 −43.77 −28.15 −29.58 224507-9 Liquid methanol Monte Carlo simulations J. Chem. Phys. 127, 224507 2007 R eu s e of AIP Pub lis h i ng c ontent is s ub j e c t to the te rms : http s ://pub lis h i ng.a i p.o r g/autho rs / ri ght sand - pe rmissi on s . D o w n l oaded to IP: 150.214.182.116 On: W ed, 19 O c t 2016 14:19:58