J. Phys. Chem. 1995,99, 12239-12248
12239
Molecular Mechanics Parameters for Fluorine-Substituted Methanes from ab Initio Quantum Calculations Bruce J. Palmer* and James L. Anchell Pacific Northwest Laboratory, Richland, Washington 99352 Received: January 17, 1995; In Final Form: June 8, 199P
Molecular mechanics parameters for the homologous series of fluorine-substituted methanes, CHxF4-x, are derived from ab initio Hamee-Fock surfaces. Ab initio intermolecular potential energy surfaces are calculated using a 6-3 1+G* basis set and include correlation using second-order Mgller-Plesset perturbation theory. A least squares fit of the ab initio surface to a standard molecular mechanics potential function, including Lennard-Jones interactions plus partial charges, is then performed. The thermodynamic properties of the resulting molecular mechanics potential are calculated using conventional molecular dynamics simulations and compared to experimental results. Additional fine-tuning of the molecular mechanics potential is then performed to optimize the agreement with experimental results. The effect of including high-energy configurations in the fit is systematically investigated.
I. Introduction The expanding availability of high quality quantum mechanical calculations for reasonably large systems has prompted increased interest in determining molecular mechanics parameters directly from points on the intra- and intermolecular potential energy surface. Most previous determinations of molecular mechanics parameters have relied heavily on comparison with experimental data, but this approach has several drawbacks. The most obvious is the lack of large amounts of experimental data for many of the model systems used in the determinations. Even when experimental information is available, it is frequently difficult to sort out how the different parameters in the potential function affect various experimental properties. This is particularly true for intermolecular forces, where there are relatively few experimental probes capable of examining individual interactions. For most systems, the quality of the intermolecular potential energy surface can only be checked by comparison with the bulk thermodynamic properties of the corresponding liquid or solid. Bulk properties, such as the enthalpy of vaporization, heat capacities, and pressuredensity-temperature behavior, usually depend in a complicated way on the potential energy surface as a whole. If there is a significant discrepancy between the model and experiment, it is difficult to determine what parameters need to be adjusted in order to achieve better agreement. Provided that they can be done to sufficient accuracy, ab initio quantum calculations offer one solution to the lack of detailed experimental information on intermolecular interactions. Points on the intermolecular potential energy surface can be obtained directly from the calculations and used to fit parameters in a molecular potential function. If necessary, these parameters can then be further fine-tuned to match experimental data. Although several recent papers have described general algorithms for extracting molecular mechanics parameters from quantum calculations,’-3 relatively few attempts to fit parameters solely from ab initio intermolecular potential energy surfaces have been reported in the literature. One of the earliest efforts to derive an intermolecular potential energy function from quantum calculations was the MCY potential of Matsuoka et aL4 This potential was obtained by fitting a large number of @
Abstract published in Advance ACS Absrracrs, July 15, 1995.
0022-365419512099-12239$09.0010
points on the water dimer surface calculated at the configuration interaction level. However, subsequent molecular simulations based on this potential showed that the MCY model gave relatively poor results for the thermodynamic properties of watere5 More recently, ab initio quantum calculations have been used to derive intermolecular potential functions for methane and nitr~gen.~.’Bohm et al. have proposed a combination of ab initio calculations and fits to experimental data to develop intermolecular potentials.8 This approach has been applied to several of the systems investigated here, although the authors do not report any results for the bulk properties of the corresponding liquids. This paper will focus on using quantum calculations to develop parameters for fluorine-substituted methanes. These compounds form a homologous series, starting with methane and continuing through tetrafluoromethane. The relatively small size of these molecules and the technological importance of several members of this series make them attractive candidates for both quantum calculations and molecular dynamics simulations. The potential functions derived for the fluorinated methane series can also serve as a basis for developing parameters suitable for simulations of other fluorine-substituted hydrocarbons. Many such molecules are of industrial and technological interest and include compounds used as refrigerants, surfactants, and wetting agents. 11. Methodology The development of the intermolecular potential energy functions begins with a quantum calculation of the minimumenergy ground-state structure for the fluorinated methane monomer. A series of quantum calculations on the dimer is then performed to map out points on the intermolecular potential energy surface. The internal geometry of the individual molecules in the dimer is assumed to remain the same as the optimized monomer geometry. The various points on the dimer potential energy surface are then used to fit the parameters in an analytic expression for the dimer potential energy. The result of the quantum calculations on the different configurations is a set of interaction energies Ea for each configuration a. Depending on the symmetry of the molecule under study, seven to ten orientations of the dimer are included in the configurations, with each orientation containing ap0 1995 American Chemical Society
Palmer and Anchell
12240 J. Phys. Chem., Vol. 99, No. 32, 1995
proximately 20 values of the carbon-carbon separation. The total number of calculated Ea per dimer range from about 150 to 200. The dimer geometry for configuration a is described by the generalized coordinate vector &. The goal is to represent the dimer interactions by a function @(R,a),where R describes the dimer geometry and a is a vector of adjustable parameters. The values of a can be found by minimizing the objective function
with respect to the individual elements in a.2%3 The weights wa can be chosen so that certain parts of the potential energy surface are favored more heavily in the fit than others. For these calculations, the function @(R,a)was assumed to have the form
The weights w a were chosen so that
w, = 1 E, < E,,, w,= 0 E, > E,,,
This particular choice ensures that no configurations with energy greater than Ecutare included in the fit. For each molecule, fits were done for several values of E,,, to provide some estimate of the sensitivity of the results on the sampling of the configurations. Once the curve fits were obtained, molecular dynamics simulations were performed at a liquid density corresponding to a point along the vapor-liquid coexistence curve. The pressure of the liquid and the enthalpy of vaporization were calculated from the simulation and compared to the experimental values. The enthalpy of vaporization is given by the formula M v a p
where i runs over all atoms in molecule 1,j runs over all atoms in molecule 2 , and ru is the interatomic distance between centers i and j . The Lennard-Jones function &(r& is given by
= 0 rij > r,
(3)
where the parameters av and bv are determined by the requirements that both 4~and its gradient vanish at the cutoff distance r,. The values of AJ and Bi are determined from the mixing rules Aij = (AiAj)1/2and Bv = (BiBj)”*. The Ai and Bi, along with the partial charges qi, are parameters associated with each of the distinct atomic species in the molecule. The Ai, Bi, and qi are obtained by direct minimization of x2. The qi are generally constrained by the requirement that the sum of qi equals the total molecular charge. Further constraints can be placed on the qi by requiring that the qi reproduce the calculated molecular dipole moment. The calculated molecular dipole moment was used in all fits instead of the experimental value in order to minimize the amount of experimental detail needed to determine potential function parameters. As discussed above, the starting point for the determination of the dimer potential function is the calculation of the interaction energy for a representative number of different dimer configurations. How large and varied the sample must be in order to be “representative” is mostly a matter of guesswork. If there are N adjustable parameters in the potential function @, then at least N configurations are needed to determine a, but beyond that there are no fixed criteria for choosing the different configurations. For this work, several different orientations of the two molecules in the dimer relative to each other were selected, and the interaction energies for different values of the carbon-carbon distance were calculated. The maximum carboncarbon distance was at 10 A, and points were generally calculated down to distances short enough so that the dimer was high on the repulsive wall. The orientations of the two molecules were chosen to highlight different atom-atom interactions. The focus was on the fluorine-fluorine, hydrogen-hydrogen, and fluorinehydrogen interactions, but an effort was made to also sample the fluorine-carbon and hydrogen-carbon interactions. For all systems, the total number of configurations was far in excess of the number of parameters in the potential.
= Hg - HI = u, PgVg- u, - PIV,
+
(4)
where U is the average potential energy per molecule, V is the volume per molecule, P is the pressure, H is the enthalpy, and the subscripts g and 1 specify the gas and liquid phases. For systems at the vapor-liquid coexistence point, P, = PI. However, because most molecular models do not have the same phase behavior as their experimental counterparts, it is worthwhile to retain the distinction between the gas and liquid pressures. Many workers simplify eq 4 by assuming that the vapor can be described by the ideal gas equatimg This implies that U, 0 and that PgV, = RT, where R is the ideal gas constant and Tis the temperature. By use of these approximations, eq 4 becomes
-
AHqap= RT
- U, - PIV,
This is the expression that was used to evaluate AHvapin this paper. The ideal gas approximation is particularly attractive because it eliminates the need to do a second simulation of the gas phase. For most systems, the experimental vapor pressure is usually in the neighborhood of 1-10 atm, and because V, >> V, the PIV term makes a negligible contribution to AHqapand is usually ignored. The enthalpy of vaporization can then be calculated using the formula
AHqap= RT - V I If this formula is used, care must be taken to verify independently that the pressure of the liquid is close to the experimental vapor pressure. Molecular models often give pressures that are hundreds of atmospheres, if not kilobars, different from the experimental values. A large value for the model pressure can substantially alter the calculated value of AHvap. Unless otherwise specified, the results from the ab initio calculations were carried out at the MP2 frozen core level of theory in a standard Pople 6-31+G* basis.l0 To test the quality of this level of theory, the 6-31+G* basis was compared with calculations on the CHzF2 monomer using several other bases. Table 1 shows the geometry and total energies obtained from the optimization of the monomer in a series of standard Pople bases, along with the experimental values of Lide.” With respect to the optimized geometry, all of the bases predicted comparable structures, Le., the standard deviations in the values of predicted bond lengths and bond angles were on the order of 0.03 A and lo,respectively. The higher level bases are also in good agreement with the experimental geometry.
J. Phys. Chem., Vol. 99, No. 32, 1995 12241
Parameters for Fluorine-Substituted Methanes TABLE 1: Optimized Geometry and Total Energies for CHzF2 for Several Basis S e w basis
RCH
RCF
LHCH LFCF
6-31G 6-31SG 6-3 1G* 6-31+G* exptb
1.071 1.072 1.078 1.077 1.091
1.383 1.391 1.338 1.342 1.358
114.1 115.2 112.4 113.4 112.1
108.6 107.9 108.6 108.3 108.2
energy -237.823 -237.835 -237.896 -237.910
a Bond lengths and angles are given in units of angstroms and degrees and total energies are reported in hartrees. From Lide."
For this study, an accurate description of the repulsive regions of the potential energy surface was desired so that the classical turning points in the molecular dynamics simulation would be correctly predicted. The weaker basis sets were insufficient in this regard, as they suffered from large basis set superposition error (BSSE)I2 in regions where the distance between the monomer units was small and the potential energy was repulsive. This is a consequence of the large energy lowering in the core atomic regions of neighboring monomers. Calculations on the CHzFz monomer indicated that, for the 6-31+G* basis set, the counterpoise correction to the BSSE was less than 1 kcavmol in regions of the potential energy surface where the repulsive energy was 10 kcal/mol or less. The size of the 6-31+G* basis was also moderate enough to carry out hundreds of calculations on various dimer configurations at a reasonable computational expense. The largest calculations involving the CF4 dimer included 190 basis functions. Because CH4 and CF4 do not possess permanent dipole moments, it was necessary to include correlation effects to describe the attractive regions on the surface arising from induced dipole interactions. For this purpose, the calculations were carried out at the relatively inexpensive second-order level of Merller-Plesset perturbation theory (MP2). Methane and difluoromethane were investigated in more detail than the other compounds. The relatively small number of electrons in the methane dimer made possible more extensive quantum calculations. The effects of increasing the level of correlation to MP4, changing the basis size, and correcting for the BSSE were also studied for this system. Some additional calculations were done for difluoromethane to estimate the effects of varying the form of the analytic potential.
111. Results The molecular dynamics calculations all consisted of short simulations, 12.5 ps in length, on 256 molecules. The simulations were performed after a 2.5 ps thermalization period. For some systems, more than one 12.5 ps interval was simulated. These additional simulations indicated that the uncertainty in the pressure was on the order of 20-30 atm and the uncertainty in the enthalpy of vaporization was on the order of 0.02-0.03 kcavmol. The temperature was controlled using a Nos6 thermostatI3 and the rigid molecular geometry was maintained using a variant of the SHAKE algorithm.I4.l5 The equations of motion were integrated using the velocity Verlet algorithm recast as a three-point predictor-correctorI6 with a time step of 2.5 fs. The Coulomb interactions were truncated, using a shifted force scheme similar to the truncation used for Lennard-Jones potential. The pairwise interactions were all of the form
= 0 rij
r,
with the cutoff distance r, set at 9.5
(6)
A.
The constants ai, and
H
n
Figure 1. Methane dimer orientations used in ab initio calculations. Points on the potential surface were calculated by varying the length of the carbon-carbon separation along the horizontal axis.
TABLE 2: Potential Function Parameters for Methane, Derived from Fits to the ab Znitio Potential Energy Surface E,", kcal/mol qc (le/) 1.0 2.0 3.0 5.0 10.0 adj
-0.416 -0.648 -0.596 -0.564 -0.820 -0.416
qH
Ail2 Bg6 Bp (lei) kca11/128, kcalt/128, kca11/68, kcal"6 8,
0.104 0.162 0.149 0.141 0.205 0.104
3.344 3.361 3.283 3.256 3.165 3.200
1.979 1.883 1.892 1.788 1.834 1.910
3.096 3.430 3.306 2.885 2.672 3.200
1.390 0.000 0.000 0.000 0.000 1.390
bo are again determined by the requirement that & and its gradient vanish at the cutoff distance r,. To check the pressure calculation, the code was run in a constant pressure mode using Andersen's extended Hamiltonian algorithm." The conservation of the associated Hamiltonian was comparable to that seen for constant temperature simulations, indicating that the pressure calculation was correct. A. Methane. The dimer orientations used for the methane calculations are shown in Figure 1. The length of the carbonhydrogen bond as calculated by optimizing the geometry of the monomer was 1.086 A, which is in good agreement with the experimental value of 1.093 A.'* The potential energy was calculated for different values of the carbon-carbon separation, and the resulting configurations and potential energy values were used in eq 1 to obtain values for the Lennard-Jones parameters and the partial charges. Fits to the ab initio surface were performed for several values of Ecutso that the behavior of the parameters as a function of the energy cutoff could be examined. The results are summarized in Table 2. The parameters A!"' and can be interpreted as the interaction distances where the contribution for the corresponding term in the Lennard-Jones potential reaches a value of 1 kcavmol. From Table 2 it can be seen that the partial charges and the dispersion coefficients Bc and BH are sensitive to the value of .Ecut,while the coefficients of the repulsion terms, Ac and AH. are relatively stable. The dispersion coefficient for the hydrogen atom, BH,disappears entirely for E,,, greater than 1.0 kcaymol. The dispersion coefficient BH also vanishes for several of the fluorinated methanes. The fact that the dispersion coefficients are oc-
Palmer and Anchell
12242 J. Phys. Chem., Vol. 99,No. 32, I995
TABLE 3: Thermodynamic Properties of the Potential Functions for Methane Listed in Table 2a E,,, kcal/mol
P atm
AH,,, kcaVmol
UE kcal/mol
Na
1.o 2.0 3.0 5.0 10.0 adj exptb
895 99 I 1123 1571 1705 6 10.3
-0.04 -0.28 -0.55 - 1.46 -1.69 1.45 1.58
0.046 0.088 0.096 0.234 0.301
96 100 101 106 1 IO
Also listed are the root mean square deviation u~ of the potential function from the ab inirio points and the total number of points Na used in the fit. From ASHRAE.'9
casionally zero is the reason that the Lennard-Jones potential is described using the coefficients A0 and Bu instead of the more conventional form
where €0 is the well depth and q is the hard-sphere radius. There is no combination of EU and a0 that will give a nonzero coefficient for the r-I2 term and zero for the r-6 term. Attempts to fit the potential energy surfaces using c and u were generally unsuccessful. The results of simulations on the different potential models in Table 2 are shown in Table 3, along with the experimental values.I9 The simulations were all done at 150 K and a density of 0.358 g/cm3,which corresponds to the density of the saturated liquid at this temperat~re.'~ Also included in Table 3 are the quantities UE and Na. The value of UE is defined as
and Nais the total number of configurations included in the fit, which increases as Ecutincreases. The value of UE is a measure of how much the fitted potential energy surface deviates on average from the ab initio values. For the pressure and AHvap, the agreement with experimental data for all the potential functions is fairly poor, but several trends can be observed. The most important of these is that the pressure increases and AHvap decreases as Ecutis increased. This behavior is characteristic of all the compounds studied in this paper. The results also show that the deviation UE increases rapidly as Ecutis increased, even though the total number of configurations increases by only a small amount. Increasing Ecutfrom 1.O to 2.0 kcaYmo1 adds only four configurations to the data set used in the fit but causes the deviation between the calculated potential surface and the fit to nearly double. This suggests that the function (2) does not represent the repulsive part of the potential surface very well. By use of the results of the curve for E,,, = 1.0 kcaYmol as a starting point, each of the parameters in the potential function was varied individually to obtain an estimate of how the pressure and enthalpy of vaporization varied as a function of the parameters. This information was then used to construct a potential function that closely matched the experimental results. The parameters for this "adjusted" function are also listed in the last row of Table 2. Within the uncertainties of the calculation, the function reproduces the pressure and is within 10% of the experimental heat of vaporization. However, it does not appear possible to improve agreement beyond this point. The reason for this can be seen in Figure 2, where the enthalpy of vaporization is plotted as a function of pressure for all the
-
4
1.6
t-
1.2
1
j
Y
0
8 4
0.8
-
0.4
-
$
0.0
0 0 0% 0
>
-200
0
200 400 Pressure (atm)
600
800
Figure 2. Enthalpy of vaporization as a function of pressure for different parameterizations of the Cb-CHd potential: (0)simulation; (a)experiment. different potential parameters investigated in this study. Many more parameter sets were investigated than listed in Table 2. Although the parameters were (adjusted independently, the enthalpies of vaporization all fall on a smooth curve. For methane, this curve passes close to the experimental point but does not appear to pass through it, implying that an attempt to exactly reproduce the enthalpy of vaporization will result in an increased error in the pressure. Because the methane dimer contains a relatively small number of electrons, it was possible to investigate the limitations of the electronic structure calculations in more detail than for the fluorine-containing compounds. Additional calculations were performed on methane using the 6-31++G** basis with correlation to the MP2 level, using the 6-31+G* basis with correlation to the MP4 level, and including the counterpoise correction (CC) to the BSSE for the 6-31+G* calculation at the MP2 level. The results of fits using an energy cutoff of 1.0 kcdmol are shown in Table 4, and the properties of methane calculated from these parameters are listed in Table 5. From Table 4, it can be seen that the coefficients of the repulsion terms are insensitive to changes in the quantum calculations. The dispersion coefficient for carbon shows some variation, while the dispersion coefficient for hydrogen shows considerable differences from one set of calculations to the next. The partial charges are also sensitive to the details of the calculation. The properties of the different sets of potential parameters, Table 5 , indicate that the counterpoise correction makes the potentials appreciably more repulsive, generating much higher pressures. This would be expected, because the effect of the BSSE is to increase the depth of potential wells. Increasing the correlation from MP2 to MP4 improves the properties slightly, as does using a bigger basis set. B. Fluoromethane. The orientations used for the quantum calculations on fluoromethane are shown in Figure 3. The optimized geometry of the monomer was found to be RCH= 1.090 A, RCF= 1.407 8,and LHCF = 108.051", which compare favorably with the experimental values of RCH= 1.095 A, RCF = 1.382 A, and LHCF = 108.047".20 The calculated dipole moment was 2.37 D, which is higher than the experimental value of 1.85 D.21 During the fits to the potential energy surface, the partial charges were constrained by the requirements that the total charge on the molecule vanish and that the charges reproduce the calculated dipole moment. The parameters derived from the curve fits to the energy surface are listed in Table 6. The partial charges on the fluorine and the coefficients of the repulsive terms in the Lennard-Jones potential vary weakly as a function of E,,,. The partial charges on the hydrogens and carbon are slightly more sensitive to E,,,, and
Parameters for Fluorine-Substituted Methanes
J. Phys. Chem., Vol. 99,No. 32, 1995 12243
TABLE 4: Comparison of the Potential Parameters for Methane Calculated for Ecut= 1.0 kcavmol, from ab Initio Surfaces Obtained at Dflerent Levels of Theory 6-31+G* MP2 6-31+G* MP4 6-31+G* MP2 CC 6-31++G** MP2
-0.416 -0.232 -0.500 -0.644
0.104 0.058
0.125 0.161
3.344 3.345 3.327 3.379
TABLE 5: Thermodynamic Properties of the Potential Functions for Methane Listed in Table 4 description P atm AH,,, kcal/mol 6-31+G* MP2 895 -0.04 6-31+G* MP4 5 19 0.63 6-3 1+G* MP2 CC 1944 - 1.79 6-31++G** MP2 404 0.91 the dispersion coefficients for all three atoms change dramatically as Ecutincreases. The results of simulations on fluoromethane using the different potential function parameters in Table 6 are listed in Table 7. The simulations were run at 203 K at the liquid coexistence density of 0.861 g / ~ m ~The . ~ pressures ~ are generally in better agreement with the experimental values23 than for methane. Unfortunately, no experimental data on the enthalpy of vaporization are available for fluoromethane, so it is impossible to determine whether or not both AHvapand the pressure can be fitted simultaneously. As was true for methane, the pressure and enthalpy of formation are closely correlated. A plot of AHvapas a function of pressure for all the different parameter sets investigated in this study is shown in Figure 4. The points all fall on a smooth curve, but it is not possible at present to determine whether the curve passes near the experimental point. The deviation DE of the fitted surface from the ab initio surface increases with increasing values of Ecut,similar to the behavior seen for methane. However, comparison of the fitted and ab initio surfaces for several of the orientations shows substantial differences between the surface corresponding to E,,, = 1.0 kcal/mol and the ab initio surface. For several of the orientations that are purely repulsive, a minimum forms in the fitted surface. An example of this is shown in Figure 5 for orientation 111. The minimum appears because not enough points are included on the repulsive surface to prevent its formation. Increasing E,,, to 2.0 kcdmol causes most of these extra minima to disappear and gives a surface that is in better overall qualitative agreement with the ab initio points, although a minimum still remains for configuration VI. Because the fit for E,,, = 2.0 kcal/mol appeared to match the potential much better over the entire energy range, it was chosen as the starting point for adjusting the potential to match experimental data. The pressure was found to be quite sensitive to the parameter Bc, and a slight adjustment of this parameter brings the model pressure in line with experiment. The effect on the potential surface of adjusting this parameter is also shown in Figure 5. The adjustment results in only a small change from the curve for E,,, = 2.0 kcal/mol. C. Difluoromethane. The orientations used for the difluoromethane dimer calculations are shown in Figure 6. The optimized geometry of the monomer was found to be RCH= 1.089 A, RCF = 1.375 A, LHCH = 114.477', and LFCF = 108.422'. The calculated dipole moment was 2.27 D, which is appreciably higher than the experimental value of 1.96 D." The parameter sets calculated using different values of E,,, are listed in Table 8. Simulations using these parameters were performed at 298 K and the liquid saturation density of 0.959 g/cm3.I9 The results of the simulations are listed in Table 9. The behavior
1.979 2.01 1 1.988 2.000
3.096 2.794 3.172 3.311
1.390 1.703 0.000 1.407
H)-H;kH . HkF ) 1 I
F
II
H
H
F
i
Figure 3. Fluoromethane dimer orientations used in ab initio calculations. Points on the potential surface were calculated by varying the length of the carbon-carbon separation along the horizontal axis. One of the fluorines in orientation VI11 is hidden by the out-of-plane hydrogen. of the pressure and AHvapparallels the behavior already seen for methane and fluoromethane. The pressure increases and the enthalpy of vaporization decreases as E,,, increases. The quality of the fit also decreases rapidly for higher values of E,,,. The value of DE,also listed in Table 9, rises rapidly with Ecut, even though only a few more configurations have been included in the fit. The parameters were adjusted to improve the value of the pressure, and the resulting adjusted parameter set is also listed in Table 8. The corresponding pressure and AHvapare included in Table 9. The value of AHvapfor the adjusted parameter set is about 0.5 kcdmol higher than that from experiment, which was the best that could be achieved with this potential function. This can also be seen in Figure 7, where the enthalpy of vaporization is plotted as a function of pressure for all the different parameter sets investigated in the course of this work. The points all fall on a smooth curve that misses the experimental point by a significant amount. Also included in Figure 7 are some points calculated using different exponents for the repulsive term in eq 6. Instead of the usual value of 12, the exponents 9 and 10 were also tried. The fits of these potential functions to the potential energy surface were similar to the behavior of the fits using the exponent 12, with the quality of the fit dropping rapidly as the value of E,,[ increased. The
Palmer and Anchell
12244 J. Phys. Chem., Vol. 99,No. 32, 1995
TABLE 6: Potential Function Parameters for Fluoromethane, Derived from Fits to the ab Initio Potential Energy Surface E,,, kcaUmol
1
qc ([el)
A:'2 kcali/I28,
qH (le[)
0
4.0
4
3.5
A;'" kcal'/128,
)
0
-
00
9
i v A y
-
H
Figure 4. Simulation (0)of enthalpy of vaporization as a function of pressure for different parameterizations of the CH3F-CH3F potential.
~
3
5
6
7
(::k
Figure 6. Difluoromethane dimer orientations used in ab initio
" " " " " " " ~ ' " '
4
"I:
F
0
3
B i 6 kcalII68,
F
8
3.0 2.5
kcalI168,
BY6 kcal'/6 8,
0
-
3.0
A:/'' kcal'lI2 8,
-
4.5 0
aE
2>
qF (lei)
8
9
1
0
R (A, Figure 5. Plot of the CH3F-CH3F potential surface as a function of the carbon-carbon separation distance for orientation 111: (0) ab initio surface; (0) fitted surface for Ecu,= 1.0 kcaymol; (A)fitted surface for E,,, = 2.0 kcaymol; (0)adjusted surface.
values of the pressure and density for these exponents lie on the same curve as the values generated by the exponent 12, indicating that the discrepancy in AHvapcannot be fixed by small, qualitative changes in the potential. To get better agreement with experiment, it may be necessary to include quantum effects or the effects of polarization. D. Trifluoromethane. The orientations used for the quantum calculations on trifluoromethane are the same as those shown in Figure 3, with the hydrogen and fluorine atoms
switched. The optimized geometry of the monomer was found to be RCH = 1.088 A, RCF= 1.351 A, and LHCF = 110.580", which compare well with the experimental values of RCH= 1.098 A, RCF= 1.332 A, and LHCF = l10.1".24 The value of the dipole moment was 2.016 D, which is substantially higher than the experimental value of 1.65 D.25 The parameters derived from the curve fits to the energy surface are listed in Table 10. The partial charges and repulsion coefficients are relatively constant as Ecutincreases, but the dispersion parameters vary widely. Both the carbon and hydrogen dispersion coefficients switch to zero for certain values of .Ecut,and the dispersion coefficient for fluorine fluctuates over a fairly broad range of values. Simulations were run at a temperature of 273 K and at the liquid saturation density of 1.0298 g/cm3;I9 the results are summarized in Table 11. Although the calculated pressures are relatively close to the experimental vapor pressure for all values of Ecut,the calculated values of A H v a p are consistently too high. A plot of AHvapas a function of the pressure, shown in Figure 8, shows that the discrepancy is quite large. The points all lie on a smooth curve that misses the experimental point by a wide margin. If the potential is adjusted to match the experimental pressure, the calculated value of AHvapis 0.74 kcal/mol too large. This is a 31% error in AHvap. The fact that the points all lie on the same curve again suggests that it is not possible to improve the potential by small adjustments in the parameters. Assuming that the ab initio surface is reasonably close to the true surface, the poor agreement between the calculated pressure and enthalpy of vaporization again suggests that important
Parameters for Fluorine-Substituted Methanes
J. Phys. Chem., Vol. 99, No. 32, 1995 12245
TABLE 8: Potential Function Parameters for Difluoromethane, Derived from Fits to the ab Initio Potential Energy Surface E,,, kcaVmol 1.o
qc (le[) q F (le])
0.054 0.006 0.290 0.154 0.652 0.050
2.0 3 .O 5.O
10.0 adj
-0.181 -0.171 -0.231 -0.202 -0.308 -0.180
qH ([el)
A:'* kcal1/I28,
A"" F kcal'/128,
A'/'' H kcal'/'*8,
Br/6 kcal'j68,
B F kcall" 8,
B p kcal'/68,
0.154 0.168 0.086 0.125 -0.018 0.155
2.998 2.778 3.306 3.426 2.919 2.900
2.657 2.669 2.357
1.958 1.937 1.893 1.714 1.781 1.950
3.363 3.028 4.023 3.755
2.389 2.513
O.Oo0 O.Oo0 O.Oo0 O.Oo0 O.Oo0 O.Oo0
2.096
1.974 2.650
TABLE 9: Thermodynamic Properties of the Potential Functions for Difluoromethane Listed in Table 8 E,,, kcal/mol P atm AHvaDkcaUmol DE kcaVmol N, 1.o 210 3.45 0.159 103 111 2.0 319 3.12 0.216 3.0 367 2.84 0.390 114 5.0 513 2.22 0.725 118 10.0 1218 0.39 1.183 122 adj 22 3.87 exptb 16.7 3.39 a Also listed are the root mean square deviation UE of the potential function from the ab initio points and the total number of points N , used in the fit. From ASHRAE.I9
58
y,
4.0 0 0
3.5
0
I
4,
O.Oo0
0.000
0.648
0.000
3.200
2.500
for tetrafluoromethane follow the same behavior seen in the other fluorinated methanes. The quality of the fit degrades rapidly as the value of E,,, is increased. The behavior of the model properties as a function of the parameters was investigated by varying the individual parameters about the values obtained using E,,, = 1.O kcdmol. For tetrafluoromethane, the properties are extremely sensitive to the value of the dispersion coefficient BF. A small adjustment of this coefficient is enough to bring the model in line with experimental values. While the Ecut= 1.0 kcdmol fit exhibits the worst match of properties, it remains true that good agreement with experimental values can be obtained by making only minor adjustments to this set of potential parameters. A plot of A H v a p as a function of pressure for the different parameter sets is shown in Figure 9. Once again, most of the points lie along a smooth curve, with the curve passing close to the experimental point but not through it. When the potential is adjusted so that the model pressure is close to the experimental vapor pressure, the corresponding value of AHvapis only 10% above the experimental value.
3.0
IV. Discussion
2.5
The results of the curve fits indicate that it is possible to do at least a reasonable job of reproducing the ab initio surface using the analytic form represented by eq 2. However, noticeable differences remain for all of the molecules examined. The parameters determined from the fits also show considerable variation, both as a function of the energy cutoff, E,,,, and from molecule to molecule. A major question in developing molecular mechanics force fields is whether or not the Lennard-Jones parameters for a given atomic species are transferable from one system to the next. A plot of the repulsion parameters for the adjusted potentials, A:'2, A:'2, and A:'2, as a function of the number of fluorine atoms is shown in Figure 10. The parameters A:'2 and A i l 2 are relatively constant as the number of fluorines increases and only vary over a range of 0.25 kcal'/'2 8, or less. The parameter A:'2 appears to increase steadily with fluorine number, except for a dip at difluoromethane. Even with the dip, the values of Ac seem to be centered around a value of 3.3-3.4 kcal'/I2 8,. The dispersion parameters exhibit larger variations in their values. A plot of @I2, B:'*, and B;" is shown in Figure 11. The dispersion parameter shows a broad oscillation centered around 3.2 kcal'I6 A, while B;I2 vanishes completely for difluoro- and trifluoromethane. The coefficient for fluorine, on the other hand, is relatively constant, with only small fluctuations about a value centered at 2.4-2.5 kcal'l6 A. It is also interesting to examine the behavior of the partial charges as a function of the number of fluorine atoms. These are plotted in Figure 12. Surprisingly, the charges on both fluorine and hydrogen are relatively constant, with the charge on the carbon atom adjusting itself to guarantee charge neutrality. For those cases where the dispersion coefficients are nonzero, the parameters A and B can be converted into the more commonly used e - g parameters used to characterize the
2.0 -500
-250
0
250
500
Pressure (atm)
Figure 7. Enthalpy of vaporization as a function of pressure for different parameterizations of the CH2F2-CH2F2 potential: (0)simulation using Lennard-Jones 6-12 potential; (0)simulation using LennardJones 6- 10 potential; (A)simulation using Lennard-Jones6-9 potential; (0)experiment. qualitative features, such as polarization or quantum effects, are being ignored in the model. E. Tetrafluoromethane. The dimer orientations used for the quantum calculations on tetrafluoromethane are the same as those shown in Figure 1, with the hydrogen and fluorine atoms switched. The carbon-fluorine bond distance in the optimized monomer structure was found to be 1.335 A, which compares well with the experimental value of 1.323 A.26 The potential function parameters derived from the curve fits are listed in Table 12. Like the compounds studied previously, the coefficients of the Lennard-Jones repulsion terms are relatively insensitive to the value of the energy cutoff, as is the fluorine dispersion coefficient BF. The partial charges generally increase in magnitude with increasing E,,,, except for a small dip at E,,, = 5.0 kcal/mol, and the carbon dispersion coefficient jumps between a value of approximately 2.8 kcal'/6 8, and zero. The value of the carbon repulsion coefficient Ac is noticeably larger than the values seen for the other fluorinated methanes. This may be due in part to the fact that the carbon atom is completely screened by the fluorine atoms and interactions between the carbon on one molecule and the fluorines on the other are masked by the fluorine-fluorine interactions. Simulations using these parameters were run at a temperature of 173 K and at the liquid saturation density of 1.465 g/cm3,I9 and the results are summarized in Table 13. The potential fits
12246 J. Phys. Chem., Vol. 99, No. 32, 1995
Palmer and Anchell
TABLE 10: Potential Function Parameters for Trifluoromethane, Derived from Fits to the ab Initio Potential Energy Surface EcUt kcaumol
qc (lei)
q F ([el)
qH(le[) A&/'*kcal1/I28,
1.o 2.0 3.0 5.0 10.0 adj
0.651 0.553 0.608 0.591 0.640 0.651
-0.233 -0.210 -0.223 -0.219 -0.230 -0.233
0.048 0.077 0.061 0.066 0.050 0.048
A i ' * kca11/12 8,
A i 1 *kcall/'* 8,
Bg6 kcal"6 8,
2.640 2.692 2.704 2.572 2.596 2.640
1.712 1.922 1.887 1.659 1.658 1.712
3.750 O.Oo0 1.713 3.390 3.203 3.590
3.361 2.970 2.937 3.207 3.136 3.367
TABLE 11: Thermodynamic Properties of the Potential Functions for Tfluoromethane Listed in Table UP E,,, kcal/mol P atm AH,,, kcal/mol UE kcal/mol Na 1.o 2.0 3.0 5.0 10.0 adj exptb
-115 - 107 -82 119 183 24 24.7
4.20 4.14 3.89 2.64 2.38 3.09 2.35
0.170 0.214 0.235 0.306 0.37 1
.
136 141 145 150 153
a Also listed are the root mean square deviation UE of the potential function from the ab initio points and the total number of points Na used in the fit. From ASHRAE.I9
1.o 2.0 3.0 5.0 10.0 adj exptb
-786 -633 -546 -326 -91 0 5.1
3.6
4.0
2.5 2.0
1
1
4 t 41 '
-200
0
0
85 88 90 91 95
2.4
0 0 0 @O
"8""~ 0
0
0
0
-100
0.173 0.252 0.276 0.292 0.582
8
I:
3.5 n
4
6.23 6.15 5.97 5.29 3.39 2.64 2.40
Also listed are the root mean square deviation UE of the potential function from the ab initio points and the total number of points Nu used in the fit. From ASHRAE.I9
4.5
3.0
0.000 1.944 1.781 0.000 0.000 0.000
2.237 2.721 2.687 2.174 2.246 2.237
TABLE 13: Thermodynamic Properties of the Potential Functions for Tetrafluoromethane Listed in Table 1 2 O EcutkcaUmo1 P atm AHvaakcaUmol UE kcaVmol N,
5.0
X'
Bi6 kcal'l6 A Bi6 kcall" A
100
2.0 -200
200
' '
'
' '
-120
Pressure (atm)
Figure 8. Enthalpy of vaporization as a function of pressure for different parameterizations of the CHF3-CHF3 potential: (0)simulation; (0)experiment.
'
'
' ' '
-40 40 Pressure (atm)
'
' '
I
' ' '
120
1
200
Figure 9. Enthalpy of vaporization as a function of pressure for different parameterizations of the CF4-CF4 potential: (0)simulation; (0)experiment.
TABLE 12: Potential Function Parameters for Tetrafluoromethane,Derived from Fits to the ab Initio Potential Energy Surface
1.0 2.0 3.0 5.0 10.0 adj
0.716 0.988 1.036 1.008 1.192 0.716
-0.179 -0.247 -0.259 -0.252 -0.298 -0.179
3.604 3.354 3.242 3.413 3.385 3.604
2.586 2.625 2.649 2.585 2.497 2.586
2.883 0.000 0.000 2.715 2.857 2.883
2.578 2.732 2.729 2.530 2.309 2.412
-
Lennard-Jones potential. This has been done for the adjusted parameter sets for all five methane compounds, and the results are listed in Table 14. Also included are the results of Gough, Deboldt, and Kollman (GDK) from a recent study to determine fluorine and hydrogen parameters for use in molecular mechanics simulation^.^^ The variation in the values of both E and u is quite large for all three atoms. For carbon, the values of uc bracket the value of uc used by GDK. The values of EC also bracket the value used by GDK, although in this case they are all either much higher or much lower. For CH3F, where the values of uc are comparable, the values of EC differ by more than a factor of 3. However, given the fact that the values used by GDK are within the range of values determined in this study, the two sets of parameters are probably compatible with each other.
1
2 1
2 3 4 Number of Fluorines Figure 10. Repulsion parameters plotted as a function of fluorine number: (0)A:'*; (A) A i 1 * ;(0) 0
Ai".
The same cannot be said for fluorine. The parameters for fluorine are all substantially different from the values determined by GDK. The values of UF in this study are all noticeably smaller than the value of GDK, and with the exception of CHF3, the values of CF are all larger. This strongly suggests that a smaller value of UF and a deeper well E F should be used for fluorine. Finally, for the two cases where parameters for hydrogen could be determined, they are both in reasonably good agreement with those of GDK. The hydro-
Parameters for Fluorine-Substituted Methanes
J. Phys. Chem., Vol. 99, No. 32, 1995 12247
1
2 3 4 Number of Fluorines Figure 11. Dispersion parameters plotted as a function of fluorine Bi6. number: (0)B P ; (A) Bk'6; (0) 0
0.8
0.4 0 -0.4
-0.8
I 0
I
I
I
1
2
3
4
Number of Fluorincs Figure 12. Partial charges plotted as a function of fluorine number: (0) q C ; (A) qF; (0)qH.
TABLE 14: Lennard-Jones Parameters Written in r - a Form for Ad.iusted Potentials" system
CH4 CH3F CHzF2 CHF, CF4 Goughetal.
uc
EC
3.200 3.910 2.628 3.158 4.505 3.82
0.250 0.0309 0.815 0.539 0.0172 0.1094
OF
2.360 2.809 3.116 2.773 3.50
EF
0.457 0.124 0.0342 0.108 0.061
OH
EH
2.625 0.0055 2.311 0.0181
2.42
0.015
" The values of E are in kcal/mol; the values of u are in angstroms.
gen parameters for CH3F are almost identical to GDK's values. For C b , the value of CH is much lower than G D K s value, but given the small well depth, this difference is probably not serious.
V. Conclusions The results of this study on the fluorinated methanes indicate that ab initio calculations can be an important tool in developing intermolecular potential functions, but they need to be used carefully. Even with today's computer resources, the calculated potential surface is likely to differ from the true adiabatic potential surface by amounts that are significant on the scale of nonbonding interactions. This creates a degree of uncertainty at the outset of a potential function determination that cannot be fixed by any subsequent manipulation of the fitting process. The analytic form used in the fit to the fluorinated methane surface reproduces the global features of the ab initio surface reasonably well, but it is not capable of reproducing it exactly. This may be the largest source of uncertainty in the fitted parameters. The sensitivity of the parameters to the value of Ecutis probably due to the inability of the r-'* term in the
Lennard-Jones potential to model the repulsive interactions high on the repulsive wall. Further problems may occur because both the partial charges and the Lennard-Jones r-6 term are trying to account for the intermediate range interactions. As Ecutincreases, the Coulomb and dispersion interactions appear to be trading off with each other, which may account for the large variations in the partial charges and the dispersion coefficients. As more repulsive interactions are included in the fit, the repulsion coefficients adjust to minimize the difference between the fitted and ab initio surfaces at high energies. The partial charges and dispersion coefficients then adjust, trying to patch up the behavior of the fit at lower energies. For most typical molecular dynamics simulations, the behavior of the potential surface at lower energies is most relevant to predicting the behavior of the system, so it seems prudent to use parameters obtained using low values of E,,,. However, as the results from fluoromethane show, this can also lead to problems. The results of the more extensive calculations on C& indicate that substantial improvements to this method will not be found by pursuing higher levels of correlation or by increasing the size of the atomic orbital bases. Other factors that may affect the results more strongly are the inclusion of three-body interactions and intemal degrees of freedom (bond stretching, bending, and torsion) in the analytic potential. These effects could be parameterized using data obtained from trimer and monomer calculations, respectively. The variations in the parameters when different numbers of configurations are included in the fit could be minimized if analytic forms could be found that do a better job of reproducing the ab initio surface. However, this would be useful only if the overall number of adjustable parameters in the potential can be kept low. A large number of parameters in the analytic form imply a high probability that two or more parameters are describing approximately the same portion of the potential energy surface, which would again result in a high sensitivity to the choice of configurations used in the fit. An interesting discussion of alternative forms for the repulsive potential has recently been presented by Stone and Tong.28 The parameters derived from the fits can be improved by subsequent adjustment to match experimental data. It is important to consider more than one piece of experimental information when doing so. As the results here indicate, a potential which successfully matches only one experimental number can be seriously deficient with respect to others. Correlation plots, such as the pressure-enthalpy plots used here, can provide important information on how well one can expect to simultaneously match different experimental numbers. Overall, the use of ab initio calculations as a starting point to develop intermolecular potential functions for molecular mechanics applications appears promising, but many issues require further examination. Among these are the effects of the quality of the ab initio surface, the sampling of the configurations, the choice of the analytic form of the potential energy fits, and the transferability of the potential function parameters. Other areas not touched on in this paper are the inclusion of multibody effects, such as polarizability, and the inclusion of trimer and higher order clusters in the curve fits. Molecular flexibility and quantum effects also merit additional study. Acknowledgment. The authors are indebted to Dave Dixon and Anthony C. Hess for many useful discussions. This work was supported by Conservation and Renewable Energy, Office of Industrial Technologies, Advanced Industrial Concepts Program,US.Department of Energy under Contract DE-AC06-
12248 J. Phys. Chem., Vol. 99, No. 32, 1995 76RLO 1830. Pacific Northwest Laboratory is operated for the U.S. Department of Energy by Battelle Memorial Institute. J.L.A. gratefully acknowledges support from the U.S.Department of Energy, Division of University and Industrial Programs, Office of Energy Research, as a postgraduate fellow at Pacific Northwest Laboratory. Portions of this work were completed on the computer resources at the National Energy Research Supercomputer Center with a grant provided by the Scientific Computing Staff, Office of Energy Research, US.Department of Energy.
References and Notes (1) Palmo, K.; Pietila, L.-0.; Krimm, S . J . Compuf. Chem. 1991, 12, 385. (2) AlemAn, C.; Canela, E. I.; Franco, R.; Orozco, M. J . Comput. Chem. 1991, 12, 664. (3) Maple, J. R.; Hwang, M.-J.; Stockfisch, T. P.; Dinur, U.; Waldman, M.; Ewig, C. S.; Hagler, A. T. J . Comput. Chem. 1994, 15, 162. (4) Matsuoka, 0.;Clementi, E.; Yoshimine, M. J . Chem. Phys. 1976, 64, 1351. (5) Owicki, J. C.; Scheraga, H. A. J . Am. Chem. SOC. 1977, 99, 7403. (6) Hay, P. J.; Pack, R. T.; Martin, R. L. J . Chem. Phys. 1984, 81, 1360. (7) Schindler, H.; Vogelsang, R.; Staemmler, V.; Siddiqi, M. A,; Svejda, P. Mol. Phys. 1993, 80, 1413. (8) Bohm, H. J.; Ahlrichs, R.; Scharf, P.; Schiffer, H. J . Chem. Phys. 1984, 81, 1389.
Palmer and Anchell (9) Jorgensen, W. L.; Swensen, C. J. J . Am. Chem. SOC. 1985, 107, 569. (10) Hehre, W. J.; Radom, L.; Schleyer, P. v. R.; Pople, J. A. Ab Initio Molecular Orbital Theory; Wiley: New York, 1986. (11) Lide, D. R., Jr. J . Am. Chem. SOC.1952, 74, 3548. (12) Boys, S. F.; Bemardi, F. Mol. Phys. 1970, 19, 553. (13) Nose, S. J . Chem. Phys. 1984, 81, 511. (14) Palmer, B. J. J . Comput. Phys. 1993, 104, 470. (15) Palmer, B. J.; Garrett, B. C. J . Chem. Phys. 1993, 98, 4047. (16) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon Press: Oxford, 1987. (17) Andersen, H. C. J . Chem. Phys. 1980, 72, 2384. (18) Gray, D. L.; Robiette, A. G.Mol. Phys. 1979, 37, 1901. (19) 1993 ASHRAE Handbook: Fundamentals; Parsons, R. A,, Ed.; American Society of Heating, Refrigerating and Air-conditioning Engineers: Atlanta, 1993. (20) Duncan, J. L. J . Mol. Struct. 1970, 6, 447. (21) Marshall, M. D.; Muenter, J. S . J . Mol. Spectrosc. 1980, 83. 279. (22) Shinsaka, K.; Gee, N.; Freeman, G. R. J . Chem. Thermodyn. 1985, 17, 1111. (23) Michels, A.; Wassenaar, T. Physica 1948, 14, 104. (24) Gosh, S . N.; Trambarulo, R.; Gordy, W. J . Chem. Phys. 1952,20, 605. (25) Sutter, H.; Cole, R. H. J . Chem. Phys. 1970, 53, 132. (26) Bowen, H. M. J. Trans. Faraday SOC.1954, 50, 444. (27) Gough, C. A.; Debolt, S . E.; Kollman, P. A. J . Comput. Chem. 1992, 13, 963. (28) Stone, A. J.; Tong, C.-S. J . Comput. Chem. 1994, 15, 1377.
JP9501658
