Intermolecular Potentials of Methane Assessed by Second Virial

Subscriber access provided by HKU Libraries

Article

Potentials of Methane Assessed by Second Virial Coefficients, Ab Initio Dimer Interaction Energies and Aggregate Cohesive Energies Douglas da Siva Ribeiro J. Phys. Chem. A, Just Accepted Manuscript • Publication Date (Web): 04 May 2017 Downloaded from http://pubs.acs.org on May 5, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry A is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Potentials of Methane Assessed by Second Virial Coefficients, Ab Initio Dimer Interaction Energies and Aggregate Cohesive Energies

Douglas S. Ribeiro Petroleo Brasileiro S.A. – RGN- Comércio de Produtos Claros – Regional São Paulo Avenida Paulista, 901 –Cerqueira Cesar – São Paulo CEP 01311-100 e-mail: [email protected]

ACS Paragon Plus Environment

1


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 31

ABSTRACT

This study presents computations of three energy related properties for 26 previously published multi-site intermolecular potentials of methane: MM2, MM3, MM2en, MM3en, MM2mc, MM3mc, MM3envir, RMK, OLPS all-atom, MUB-2, AMBER, BOYD, Williams, Sheikh, MG, Tsuzuki, E2-Gay, E4-Gay, MP4exp-6(iii), MP4exp6(iv), Rowley-A, Rowley-B, TraPPE-EH, Ouyang, CLC and Chao and three united atom potentials: Saager-Fischer (SF), OPLS united atom and HFD. The three properties analyzed are the second virial coefficients for 14 temperature points in the range of 110 to 623.15 K, the interaction energies for 12 orientations of the methane dimer as a function of distance followed by a comparison to three ab initio data sets and the cohesive energy of the aggregate of 512 methane molecules. The latter computed energies are correlated to latent heat of evaporation of 11 potentials and are proposed as surrogate approximate parameters for ∆Hvap for the studied potentials. Ten best performing potentials are selected by rms order in each one of the properties and three of them are found to be present simultaneously in the three sets: Tsuzuki, MM3mc and MM2mc. Based on the cohesive energy of the aggregate a quantitative measure of the anisotropy of the potentials is proposed. The results are discussed on the basis of anisotropy, non-additivity and ability of the potentials to reproduce ab initio data. It is concluded that the non-additivity of the pair potentials holds and the available ab initio data did not lead to pair potentials that are cohesive enough to reproduce accurately the second virial coefficients.


2

Page 3 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


INTRODUCTION

The intermolecular interaction potential of methane has been studied for decades now. The potential of Williams1, which has received considerable attention, dates back to 1967. This continued interest has been ascribed to the importance of non-polar interactions in biological systems. Methane is the first homologue of the all-important hydrocarbon class and can therefore be an excellent model to study the dispersion energy in higher alkanes. It is also the principal component of natural gas with the consequent interest of the oil industry2,3 in developing good simulation tools to study its mixture with carbon dioxide and water. Perhaps one of the main reasons why so many studies have been devoted to methane is related to its small size. Deep, intricate theoretical computational studies become feasible with small molecules with the added advantage that there is a huge amount of thermophysical data available4. The intermolecular potentials are needed in computer simulations by both molecular dynamics (MD)5 and Monte Carlo (MC)6 methods. Early potentials, so called united atom models because the CH4 molecule is represented by a single atom, have attained quite successful results despite being so simple. The Saager-Fischer (SF)7 model has been shown to reproduce quite well the pressure and heat capacities of methane at several temperatures of both liquid and gas phases. The OLPS united-atom potential8, which is quite similar to the SF model and identical to the TraPPE potential, has also been demonstrated to be successful to reproduce densities, specific heat capacities, compressibilities and pressures to a wide variety of state points2.


3


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 31

Despite its simplicity which favors computer simulations, these united-atom potentials are isotropic and cannot discern the energy of different orientations of methane dimers as obtained by ab initio methods. A rather strong anisotropy of the potential has been acknowledged by many groups9,10,11. Accordingly, some authors12,13 have pointed out that an all atom model, four hydrogens plus carbon, should describe better the physicochemical behavior of methane. Potentials obtained by fitting a parametric function to bulk liquid or solid properties are known as effective potentials as they embody the true pair potential plus the many body effects present in the condensed phase14. Crystal sublimation energies1,15, Latent heat of vaporization and liquid densities7,8,12, pressure second virial coefficients16 and nonbonded interactions in conformational analysis of hydrocarbons17,18 are all examples of properties used to derive empirical intermolecular potentials. Pair potentials are obtained by ab initio methods. Large, flexible basis sets have to be used along with Moller Plesset perturbative methods (MP2, MP3, MP4) to calculate the dispersive interaction between dimers. Many research groups9,10,19,20,21 have published pair potentials from ab initio data and applied these potentials to computer simulations to ascertain its adequacy. It is normally taken for granted that the potentials are pair additive and can be used directly in condensed phases. It is a frequent feature of the studies in this field to run MD or MC simulations both with the proposed potential and literature potentials to compare them. Self-diffusion coefficients, velocity auto correlation functions and radial distribution functions are common properties used to test potentials and compare them.


4

Page 5 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The functional form of the potential and the strategy to parameterize it are also the subject of various different approaches among authors. Most authors adopt either the 12-6 Lennard Jones, the Buckingham function or, less frequently, the Morse potential. The strategy to parameterization is an art on its own. The number of independent parameters, the use of previously determined parameters, the use of weights to calculate rms deviations, determining parameters stepwise following some prerequisites or all at once, the number of properties used in the parameterization, the choice of temperature and density when bulk properties are used, the commitment or not to transferability to larger molecules are all variables that have been dealt with by a number of authors8,9,10,13,19,21,22. Nagy et al.23 published a comprehensive study of ten empirical force fields analyzing, inter alia, pressure second virial coefficients, dimers energies, latent heats of vaporization (∆Hvap) and liquid densities at the boiling point of methane by Monte Carlo simulations. The authors also presented ten novel potentials parameterized using combinations of these physical properties. They noted a clear trend within the studied potentials: except for MM2, all of them underestimated the second virial coefficients while yielding reasonable ∆Hvap values. The authors ruled out non-additivity as the explanation of this finding. They relied on the work of Szczesniak et al.11 on the ab initio calculation of methane trimers which concluded that the three body potential is very small. Given the diversity of different potentials, functional forms, fitting strategies, limited number of analyzed potentials in each simulation study, properties fitted and thermodynamic states involved there is a need to compare these potentials using a common method and criterion.


5


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 31

The main objective of this study is to present such a means of comparing the most known potentials by analyzing energy related properties as the second virial coefficients (B(T)) in the temperature range of 110 – 623 K, dimers energy and aggregate cohesive energies. Three united-atom potentials: Saager-Fischer (SF)7, OPLS8 and HFD24 and 26 multi-site potentials: from Chao et al.9 CLC, from Nagy et.al.23: MM2en, MM3en, MM2mc, MM3mc and MM3envir,. from Rowley and Pakkanen22: Rowley-A, Rowley-B, RMK16, OLPS all-atom13, Ouyang20, from Hayes et.al.10 : MP4exp-6(iii), MP4exp-6(iv), Chao25, Tsuzuki21, MM217, MM315, from Gay et al.19 E2-Gay and E4-Gay, Boyd26, AMBER27, MUB-218, Williams1, Sheikh28, MG29, and TraPPE-EH12-were studied. These potentials have been chosen on the basis of frequency of citations in the literature related solely to methane or because they have been published recently. The calculations of the second virial coefficients presented by Nagy et al.23 are expanded in this study to all the mentioned potentials in a wider temperature range. The energy of methane dimers in twelve different orientations are calculated by each studied multi- site potential as a function of distance and the results compared to three literature ab initio computations9,10,21. The attractive part of the resulting curve is fitted to a Lennard Jones function so that the potential well depth and the position of the minimum are obtained. The Saager-Fischer (SF) potential is used in a liquid simulation at the normal boiling point of methane from which points of the phase space are extracted and the carbon positions are used to construct an aggregate of 512 molecules. The energies of the aggregates are calculated by each one of the studied potentials and correlated to latent heat of evaporation obtained through MC simulations. Approximate surrogate ∆Hvap


6

Page 7 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


values are proposed. A quantitative measurement of the anisotropy of the potentials is also proposed. The top ten best performing potentials in each property out of 26 are selected based on rms deviations and it turns out that three potentials are present simultaneously at the three sets. The results of the three energy related properties are presented and discussed on the basis of the similarities among the studied potentials, their anisotropies and nonadditivity. METHODS Second Virial Coefficients. The second virial coefficients have been obtained by the numerical integration of equation 1 for the isotropic united-atom potentials: The OLPS united-atom, the Saager-Fischer and the HFD potential.

BT = −2πN 1 − exp −

Ur r dr kT 1

Na is Avogadro’s number, r is the distance between two methane molecules, U(r) is the intermolecular potential, Kb is Boltzmann constant and T is the absolute temperature. B(T) values for the multi-site model potentials have to be acquired in a way that averages all angular orientations between the two molecules. This is accomplished19,30 through equation 2 which was integrated numerically for all the multi-site potentials of this study. BT = −

! ! N R dR dθ sinθ dφ 16π


7


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

!

!

Page 8 of 31

!

sinθ dθ dφ dω $e%&'(;* − 1+

2 Molecule 1 is centered at the origin and molecule 2 is located at the polar coordinate (R, θ, φ) and rotated along its Eulerian angles θ2, ϕ2 and ω2. β = 1/kbT, α is the orientation of molecule 1 and 2. The integration was performed stepwise in 4 parts: from zero to 2.65 Å a hard-core !

potential ( , N 2.65, was assumed which contributes 23.47 cm3 / mol to the second virial coefficient. (B(T)). The second part (from 2.65 to 10.0 Å) was calculated using angular differentials of π/12 (equivalent to 15 degrees rotations) to all angular variables and radial differential of 0.05 Å. The convergence of the integral was tested up to a differential of π/18 for the path between 3.5 to 5.0 Å which resulted in a difference of only 0.6 % from the result with angular differential of π/12, thus the latter was found to be suitable. The third part (from 10.0 to 23 Å) used angular differentials of π / 6 and radial differential of 0.05 Å. The fourth part (from 23 to 40 Å) used angular differentials of π / 6 and radial differential of 0.10 Å. The long range correction19 was ignored as its contribution is less than 0.1 units. For the MM2en, MM3en, MM2mc, MM3mc, Tsuzuki, MM2, MM3, Boyd and Williams potentials the first part of the integration had to be extended to 3.0 Å contributing 34.06 cm3 / mol to B(T). For Rowley-A and Rowley-B the first part had to be extended to 3.5 Å which contributes 54.08 cm3 / mol to B(T). These functions are not defined to work at closer range. Ab Initio Dimer Energies. The interaction energies of the dimers of methane for all the studied multi-site potentials were calculated as a function of carbon-carbon distances


8

Page 9 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


from 3.0 to 8.0 Å in steps of 0.2 Å for the twelve different orientations shown in Figure 1. The obtained results were compared to the ab initio energies of Chao et. al.9, Tsuzuki et al.21 and Hayes et. al.10 and the rms values calculated through equation 4. The fitting of the energies obtained by each combination of potential and dimer orientation to a Lennard Jones function was carried out by the minimization of the residual sum of squares (RSS) as follows: The points of the attractive part of the curve plus the first positive energy point were read and the point of minimum energy found. The latter was used to provide the reference values of Sigma and Epsilon. The value of Sigma was varied with steps of 0.01 Å within the range ± 0.60 Å from the reference value and Epsilon was varied in steps of 0.002 kcal/mol within the range ± 0.20 kcal/mol from the corresponding reference values. The energy of each point in this grid was evaluated as well as the sum of the squared residuals (RSS). Whenever the RSS was lower than the previous stored value the corresponding Sigma and Epsilon were stored. The depth of the potential well (ε / kcal/mol) is the obtained Epsilon and the location of the energy minimum (r* / Å) were obtained by multiplying Sigma by 1.1224. Molecular Dynamics Simulation. The molecular aggregates used to calculate and compare the various cohesive energies of the studied potentials were obtained by extracting the carbon positions during the final production run of a molecular dynamics simulation using the potential of Saager and Fischer7. The SF potential is a united-atom model known to produce accurate thermodynamical properties in a variety of situations.5,6,7 The molecular dynamics simulation (MD) was carried out in the NVT ensemble for the normal boiling point of methane.


9


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 31

Figure 1: The orientations of Methane Dimers


10

Page 11 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


A cubic box with 512 molecules (pseudoatoms) initially arranged in a fcc structure was submitted to an equilibration period of 120.000 time steps. The production run was of 160.000 time steps. The time step was 1 fs. A cut off radius of 14.93 Å (4 σ) was used. Standard Conditions32 were applied such as: periodic boundary conditions, a fourth order predictor-corrector integration scheme, the rescaling of velocities to keep temperature constant. No long range corrections were applied to the potential energy. The runs were divided into blocks of 20.000 steps and the uncertainty quoted refers to the statistics standard deviation (1σ). The pseudoatom positions and velocities were saved at every 20.000 steps and the last 3 files were used in the production of the molecular aggregates by the addition of the hydrogen atoms. The heat of evaporation was calculated by the following equation:

∆Hvap = - Epot(liquid) + RT

3

The obtained heat of evaporation was 1.954 Kcal/mol, which exactly matches the experimental value. The pressure was -0.1 ± 9 bar. The obtained properties are in good agreement with experimental values. The heat of evaporation is of primary interest as it denotes the ability of the SF united-atom potential in modeling the behavior of methane. The objective of this simulation is just to extract the aggregates at the boiling point to form a common basis of comparison of the studied potentials.

Cohesive Energy of the Aggregates. The carbon positions extracted from the molecular dynamics simulations were used to create the aggregates of the methane molecules by adding the hydrogen atoms from a reference molecule with tetrahedral


11


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 31

symmetry and C-H bond length of 1.09 Å. The cohesive energy was optimized by the following procedure for each one of the studied potentials. A list of neighbors was created for each molecule considering a cut off radius of 13 Å. All 512 molecules in turn had their interaction energy with neighbors calculated. This energy was labeled as Initial Energy, its position was saved. That molecule was submitted to a random rotation along its center of mass, and its energy was recalculated and labeled as Final Energy. If the Final Energy was lower than the Initial Energy the loop continued to the next molecule, if it was higher the molecule returned to its saved position before the loop continued to the next molecule. After 50 rounds of such loops the energies of the entire aggregate were calculated for all of the potentials and saved. The whole process was repeated with the next potential until all of them had been used. The optimized aggregate from the previous run was used. Aggregates with modified C-H bond lengths were created as certain potentials use the interaction center shifted from the hydrogen nuclei towards the carbon nucleus. The final result of this procedure is having an aggregate being optimized by a given potential and its cohesive energy being calculated by every other potential.

Anisotropy of the Potentials. It was also of interest to calculate the cohesive energies of the aggregates in a completely randomized configuration, where the molecules have undergone random rotations: All 512 molecules were rotated randomly through their center of mass and the energy of the aggregate was calculated. This process was repeated 25 times and the results were averaged. This average energy has been termed Erand. The proposed scale of anisotropy is (Eopt – Erand)/Eopt where Eopt is the energy of the optimized aggregate. RESULTS


12

Page 13 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Second Virial Coefficients. Table 1 presents the calculated virials for the studied potentials along with the experimental data and the rms for 12 temperature points in the range of 110 – 623.29 Kelvin. MM2, MM3, Boyd, Amber, MUB-2, OLPS un-atom, OLPS all-atom, MM2en, MM3en, MM2mc, MM3mc and MM3envir are in close agreement with the corresponding values reported by Nagy et al.23 On the other hand, the rms reported by the authors could not be reproduced. In this work the rms takes its known form: 1

rms = 02 ∑2 1 $X 5 − X 678 +

4

The agreement to the reported values of E4-Gay is satisfactory considering that the authors19 used a quantum correction not used here. Even so, it is still the second best potential among those studied here based on its rms value. The results for E2-Gay have a huge discrepancy to the reported19 values and cannot be accounted for. The obtained B(T) values for the RMK potential agree closely well with the original reference16 as much as these can be checked. Unfortunately the original data have not been given in a tabular form. The united-atom potentials reproduce better the experimental virial coefficients than the multi-site potentials do and this is reflected in the overall ranking positions of 4th, HFD, 6th Saager-Fischer and 9th OLPS un-atom out of a total of 29 potentials.


13


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48

Page 14 of 31

Table 1: Pressure Second Virial Coefficientsa of Methane as a Function of Temperature Calculated by the Studied Potentials. POTENTIAL CLC MM2en MM3en MM2mcb MM3mc Rowley-A Rowley-B RMK OLPS all-atom Ouyang MP4exp-6(iii) MP4exp-6(iv) Chao Tsuzuki MM2 MM3 E2-Gay E4-Gay Boyd Amber MUB-2 Williams Sheikh MG TraPPE-EH MM3envir OLPS un-atom SF HFD EXPERIMENTALc

110 -125.1 -183.5 -187.2 -268.6 -273.2 -139.2 -259.1 -328.1 -287.9 -182.0 -284.9 -265.0 -264.5 -280.3 -693.6 -200.6 -184.2 -326.8 -247.8 -247.2 -290.6 -259.0 -248.7 -257.1 -276.1 -304.8 -277.1 -284.3 -298.9 -330

120 -103.1 -156.2 -159.5 -228.4 -232.8 -115.7 -215.5 -278.1 -242.6 -153.8 -242.7 -225.6 -222.7 -239.6 -565.2 -170.6 -155.9 -274.6 -210.7 -209.0 -247.3 -220.1 -209.7 -216.6 -235.5 -259.6 -238.3 -244.5 -255.8 -273

130 -85.6 -134.4 -137.4 -196.9 -201.2 -97.1 -182.1 -239.4 -207.6 -131.4 -209.6 -194.8 -190.2 -207.6 -472.8 -146.8 -133.4 -234.7 -181.7 -179.1 -213.5 -189.6 -179.3 -185.1 -203.6 -224.3 -207.5 -212.9 -221.9 -235

150 -59.6 -101.9 -104.4 -150.9 -154.9 -69.4 -134.5 -183.8 -157.2 -98.3 -161.4 -149.6 -143.1 -160.6 -349.7 -111.4 -100.1 -177.6 -139.0 -135.5 -164.2 -145.2 -134.8 -139.2 -156.8 -172.9 -161.9 -166.1 -172.1 -182

180 -33.8 -69.5 -71.7 -106.5 -110.0 -42.2 -89.8 -131.0 -109.4 -65.6 -114.8 -105.9 -98.0 -114.8 -243.4 -76.6 -67.2 -124.1 -97.6 -93.5 -116.7 -102.4 -91.9 -95.1 -111.2 -123.4 -117 -120.1 -123.7 -129

TEMPERATURE (KELVIN) 191.06 200 240 273.15 -26.8 -21.8 -4.8 4.9 -60.7 -54.4 -33.0 -20.8 -62.7 -56.4 -34.8 -22.4 -94.5 -86.1 -57.7 -41.7 -97.9 -89.4 -60.7 -44.5 -34.8 -29.5 -11.8 -1.7 -78.0 -69.8 -42.4 -27.2 -117.0 -107.2 -74.3 -56.0 -96.7 -87.8 -57.9 -41.1 -56.7 -50.4 -29.1 -16.9 -102.3 -93.5 -63.9 -47.2 -94.1 -85.8 -57.8 -42.0 -86.0 -77.5 -49.1 -33.0 -102.5 -93.8 -64.3 -47.6 -216.9 -198.7 -140.0 -108.6 -67.1 -60.4 -37.6 -24.7 -58.2 -51.9 -30.3 -18.1 -110.0 -100.1 -67.1 -48.8 -86.4 -78.5 -52.0 -36.9 -82.2 -74.3 -47.5 -32.3 -104.0 -95.1 -64.9 -47.9 -90.9 -82.8 -55.5 -40.1 -80.4 -72.2 -44.8 -29.3 -83.3 -75.0 -46.9 -31.1 -99.0 -90.3 -61.0 -44.4 -110.1 -100.8 -69.4 -51.6 -104.8 -96.1 -66.7 -49.9 -107.7 -98.8 -68.8 -51.7 -110.7 -101.5 -70.5 -52.9 -116.3 -105 -72 -53.3

RMS 298.15 10.6 -13.7 -15.2 -32.4 -35.0 4.1 -18.4 -45.4 -31.5 -9.8 -37.5 -32.8 -23.7 -37.9 -91.0 -17.1 -10.9 -38.2 -28.2 -23.6 -38.0 -31.2 -20.2 -22.0 -34.7 -41.4 -40.1 -41.7 -42.7 -42.82

348.15 19.1 -2.9 -4.3 -18.4 -20.9 13.0 -5.3 -29.5 -17.0 1.0 -22.9 -19.0 -9.8 -23.3 -65.2 -5.8 -0.1 -22.3 -15.1 -10.4 -23.2 -17.8 -6.7 -8.2 -20.2 -26.0 -25.4 -26.6 -27.3 -27.1

423.18 27.6 7.8 6.5 -4.6 -6.8 21.8 7.6 -13.9 -2.6 11.7 -8.5 -5.3 4.0 -8.8 -40.4 5.5 10.7 -6.7 -2.0 2.8 -8.6 -4.5 6.8 5.5 -5.6 -10.7 -10.5 -11.5 -12 -11.4

523.24 34.7 16.6 15.4 6.9 4.8 29.2 18.3 -1.0 9.4 20.6 3.5 6.0 15.5 3.2 -20.2 14.9 19.7 6.4 8.8 13.7 3.6 6.5 18.0 16.9 6.5 2.0 1.8 1.2 0.7 1.49

623.29 39.1 22.2 21.1 14.2 12.3 33.9 25.2 7.2 17.0 26.3 11.1 13.2 22.9 10.8 -7.6 20.8 25.4 14.7 15.6 20.7 11.3 13.5 25.2 24.2 14.3 10.1 9.8 9.3 8.8 9.66

103.0 68.4 66.3 26.9 23.9 94.6 38.9 2.7 20.5 71.0 18.3 28.2 32.9 19.7 157.5 60.3 69.6 4.5 37.0 39.4 15.9 31.5 40.3 36.4 22.6 8.9 19.9 16.6 10.7

a: in cm3/mol b: in bold the ten best performing multi-site potentials c : from reference 19, 31

14 ACS Paragon Plus Environment

Page 15 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The top 10 multi-site potentials based solely on the rms values are: RMK, E4-Gay, MM3envir, MUB-2, MP4exp-6(iii), Tsuzuki, OLPS all-atom, TraPPE-EH, MM3mc and MM2mc. RMK and MM3envir were modeled taking the B(T) values in the process of parametrization. MUB-2 was constructed by conformational analysis methods, OLPS all-atom was obtained by bulk liquid properties and free energy of hydration, TraPPEEH was obtained modelling bulk liquid properties and E4-Gay, MP4exp-6(iii) and Tsuzuki were obtained by ab initio dimers energy. MM3mc and MM2mc were modeled by simultaneous modeling ab initio and bulk properties. Most potentials underestimate B(T) values (less cohesive interaction) compared to the experimental points. The exceptions are the known MM2 and to a slight extent the RMK potential.

Ab Initio Dimer Eenergies. Table 2 presents the location of the energy minimum (r*) in Å and Table 3 presents the potential well depths (ε in kcal/mol) of the orientations of the methane dimers for the multi-site potentials studied.

The calculated rms presented in Table 5 for the ab initio data of Chao, Tsuzuki and Hayes were used to pick up the ten best performing potentials: MP4exp-6(iv), Tsuzuki, MP4exp-6(iii), MM3mc, MM3en, MM3envir, Williams, MM2mc, MM3 and MM2en. Except for Williams and MM3 potentials they all were obtained taking ab initio dimers energies in the process of parametrization. Nagy et. al.23 had already recognized MM3 potential as a good performer in reproducing ab initio dimer energies.


15


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 31

An inspection of Table 2 shows that the ab initio data of Tsuzuki, Chao, and Hayes agree very well with each other. It can also be seen that the potentials that best reproduce the location of the energy minimum (r*) are: MP4exp-6(iii), Tsuzuki, MP4exp-6(iv), MM3mc, MM3en, MM2mc, MUB-2, Williams, MM3 and MM3envir, basically the same that yielded the lowest rms values. The order of increasing r*, (J, I, G, H, L, K, E, F, C., D, A, B), for the twelve orientations is a common feature to most potentials including the ab initio data of Chao, Tsuzuki and Hayes. The potentials TraPPE-EH and Rowley-B show slight different orders but the real exceptions of this order are the RMK potential (L, G, H, E, F, J, I, K, C, D, A, B) and the E4-Gay potential (C, D, E, F, K, G, H, I, J, L, A, B). These two potentials disagree not only to the ab initio data and to most other potentials but also to each other. Table 2: Location of the Energy Minimum (r*) in Å of the Orientations of the Methane Dimers as a Function of the Potential. POTENTIAL A/B C CLC 5.37 4.71 MM2en 5.10 4.51 MM3en 5.06 4.49 MM2mc 4.92 4.41 MM3mc 4.86 4.39 Rowley-A 5.14 4.52 Rowley-B 4.95 4.40 RMK 4.82 4.16 OLPS all-atom 5.06 4.38 Ouyang 5.19 4.46 MP4exp-6(iii) 4.86 4.41 MP4exp-6(iv) 4.95 4.43 Chao 5.28 4.53 Tsuzuki 4.92 4.46 MM2 4.94 4.33 MM3 5.06 4.50 E2-Gay 5.15 4.40 E4-Gay 4.69 3.95 Boyd 5.12 4.47 Amber 5.13 4.46 MUB-2 4.80 4.41 Williams 4.89 4.38

D 4.75 4.53 4.51 4.43 4.40 4.55 4.42 4.19 4.42 4.49 4.43 4.46 4.59 4.47 4.37 4.52 4.42 3.98 4.49 4.49 4.42 4.40

ORIENTATION E/F G H I 4.35 3.86 3.91 3.64 4.28 3.94 3.96 3.91 4.27 4.01 4.03 3.89 4.27 3.87 3.88 3.81 4.24 3.87 3.88 3.81 4.31 4.07 4.08 3.95 4.30 3.78 3.80 3.60 4.00 3.95 3.95 4.05 4.28 3.74 3.77 3.61 4.31 3.86 3.88 3.77 4.22 3.86 3.87 3.76 4.21 3.85 3.86 3.77 4.20 3.85 3.87 3.75 4.27 3.97 3.99 3.81 4.12 3.69 3.70 3.52 4.30 3.89 3.91 3.80 4.19 4.10 4.11 4.07 4.04 4.11 4.11 4.11 4.24 3.83 3.84 3.73 4.32 3.79 3.83 3.66 4.27 3.89 3.91 3.75 4.22 3.82 3.83 3.66

J 3.48 3.87 3.86 3.78 3.79 3.91 3.64 4.03 3.59 3.73 3.73 3.73 3.70 3.78 3.48 3.76 4.06 4.11 3.68 3.64 3.74 3.65


K L 4.29 4.02 4.18 4.04 4.15 4.13 4.06 3.98 4.05 3.97 4.20 4.16 3.96 3.91 4.07 3.94 3.93 3.85 4.09 3.93 4.07 3.98 4.09 3.98 4.14 3.94 4.12 4.12 3.94 3.82 4.13 4.01 4.19 4.12 4.05 4.11 4.07 3.93 4.01 3.92 4.07 4.02 4.01 3.92

16

Page 17 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Sheikh 5.31 4.62 4.66 4.40 3.92 MG 5.29 4.61 4.65 4.39 3.86 TraPPE-EH 4.37 4.30 4.31 4.23 4.06 MM3envir 4.86 4.40 4.41 4.25 3.96 Ab initioa 4.95 4.40 4.40 4.15 3.88 Ab initiob 4.90 4.39 4.39 4.13 3.86 Ab initioc 4.98 N.A. 4.40 4.15 3.88 a: from ref. 21, b: from ref. 9, c: from ref. 10

3.95 3.89 4.06 3.91 3.88 3.87 4.02

3.77 3.69 3.97 3.83 3.75 3.73 3.73

3.70 3.64 3.97 3.81 3.68 3.67 3.67

4.22 4.02 4.16 4.00 4.15 4.15 4.06 3.98 4.19 4.02 4.19 4.01 4.20 N.A.

That the RMK potential diverges from ab initio data is understandable as it has not been developed from them but to satisfy solid state properties, but curiously the E4-Gay potential has been developed from ab initio data and yet presents a comparatively huge rms value. It is interesting to see that though the MM3en and MM2en potentials were constructed to reproduce ab initio values, they yielded shallower potential wells for all orientations. The Szczesniak11 data used by Nagy are shallower than those of later studies. The authors had already speculated that the ab initio potential well of Szczesniak were too shallow as the energy related properties were underestimated. The potentials Rowley-A, MM3, E2-Gay and Ouyang also show predominantly shallower wells than the ab initio data.

Table 3: Potential Well Depth (ε) in kcal/mol of the Orientations of the Methane Dimers as a Function of the Potential. POTENTIAL

CLC MM2en MM3en MM2mc MM3mc Rowley-A Rowley-B RMK OLPS all-atom Ouyang MP4exp-6(iii) MP4exp-6(iv)

A/B

C

D

E/F

0.428 0.140 0.141 0.187 0.190 0.119 0.204 0.141 0.227 0.170 0.171 0.158

0.263 0.205 0.208 0.273 0.272 0.193 0.284 0.332 0.298 0.245 0.256 0.243

0.291 0.201 0.205 0.270 0.271 0.190 0.284 0.324 0.295 0.246 0.250 0.238

0.279 0.253 0.249 0.300 0.294 0.226 0.303 0.428 0.337 0.306 0.306 0.296

ORIENTATION G H I 0.351 0.324 0.320 0.411 0.399 0.307 0.486 0.404 0.477 0.341 0.436 0.429

0.344 0.312 0.314 0.406 0.394 0.305 0.476 0.406 0.465 0.339 0.430 0.423

0.427 0.347 0.347 0.443 0.426 0.323 0.575 0.406 0.526 0.331 0.501 0.485


J

K

L

0.486 0.355 0.347 0.449 0.429 0.338 0.563 0.413 0.534 0.354 0.513 0.500

0.285 0.268 0.275 0.361 0.356 0.259 0.418 0.348 0.408 0.271 0.356 0.339

0.333 0.290 0.298 0.371 0.365 0.290 0.440 0.429 0.450 0.367 0.375 0.360

17


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Chao Tsuzuki MM2 MM3 E2-Gay E4-Gay Boyd Amber MUB-2 Williams Sheikh MG TraPPE-EH MM3envir Ab initioa Ab initiob Ab initioc

0.361 0.169 0.273 0.152 0.111 0.203 0.153 0.213 0.180 0.164 0.314 0.331 0.244 0.204 0.159 0.124 0.160

0.315 0.241 0.423 0.222 0.214 0.404 0.237 0.268 0.263 0.254 0.304 0.299 0.285 0.286 0.276 0.257 N.A.

0.319 0.242 0.415 0.219 0.212 0.390 0.235 0.269 0.262 0.250 0.316 0.314 0.283 0.287 0.276 0.257 0.259

0.370 0.273 0.510 0.252 0.301 0.403 0.283 0.319 0.303 0.291 0.383 0.356 0.316 0.312 0.348 0.339 0.333

0.446 0.395 0.754 0.358 0.374 0.453 0.423 0.436 0.433 0.433 0.425 0.461 0.369 0.408 0.370 0.361 0.36

0.441 0.387 0.739 0.353 0.363 0.454 0.415 0.427 0.427 0.426 0.417 0.451 0.377 0.408 0.370 0.359 0.324

0.477 0.449 0.880 0.402 0.498 0.602 0.468 0.491 0.487 0.493 0.449 0.536 0.408 0.445 0.302 0.401 0.396

Page 18 of 31

0.477 0.461 0.919 0.412 0.507 0.599 0.485 0.495 0.491 0.498 0.481 0.558 0.404 0.448 0.438 0.443 0.442

0.371 0.334 0.594 0.297 0.281 0.382 0.336 0.365 0.364 0.359 0.332 0.367 0.332 0.376 0.269 0.25 0.253

0.436 0.349 0.677 0.323 0.368 0.452 0.376 0.405 0.375 0.390 0.435 0.429 0.368 0.386 0.333 0.332 N.A.

a: from ref. 21, b: from ref. 9, c: from ref. 10 The ab initio potential wells for orientations E and F are mostly deeper than those of the studied potentials. There are 5 exceptions: RMK, E4-Gay, MG, Sheikh and MM2. For orientations G, H, I, J, K and L most of the studied potentials show deeper potential wells than the reference ab initio data.

Cohesive Energy of the Aggregates. Table 4 presents the quotients of the cohesive energies of the aggregates by the cohesive energy calculated by the reference SF potential. The potential used in the optimization is indicated at each column heading while each row heading shows the potential used to calculate the cohesive energy of that aggregate. So, for instance, the number shown at the entry whose column is labeled MM3 and row named Tsuzuki refers to the cohesive energy calculated with the potential Tsuzuki in an aggregate optimized with the MM3 potential. The energies vary only slightly along any given row. This behavior suggests that all potentials studied lead to the same optimized structure of the aggregates. All of them agree on what an optimum structure is, with two notable exceptions: When the aggregate is optimized by the E4-Gay and the RMK potential all others have their energy increased, including the E2-Gay potential.


18

Page 19 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Most potentials show greater interaction energies than the reference potential, but there are some exceptions like, for instance, CLC, MM2en and MM3en. The cohesive energies of these aggregates should indeed be greater than the reference force field if they are to reproduce well the energetic of the system. During the simulation the aggregates obviously tend to be in the optimized configuration. But due to the randomizing effect of the kinetic energy these configurations will not always be fully optimized. Thus, the potentials that show lower energy values compared to the reference are not likely to represent well the cohesiveness of the system. Table 5 presents the averages of the energy ratios contained in Table 4. The E4-Gay and RMK potentials were left out of the calculation of the averages for the reason outlined above. The small standard deviations indicate how well all these potentials agree with each other as to what is an optimized structure. As equation 3 shows, the latent heat of evaporation (∆Hvap) is directly proportional to the potential energy of the liquid phase. That is an acceptable simplification because the potential energy of the gas phase is far too small compared to that of the liquid and can be neglected. It is a known fact of liquid simulations that the potential energy of the system varies only slightly during the production run.7,32 For instance, in this study the standard deviation of the potential energy was just 0.27%. Consequently the averages of the ratios of the cohesive energy of the aggregates shown at Table 5 are tentatively used henceforth as valid surrogate values of ∆Hvap to compare the studied potentials. To do


19


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 31

that, it must be shown that these parameters do have a good linear correlation with the heat of evaporation calculated through the liquid simulations reported in the literature. The work of Hayes et al.10 presents results of ∆Hvap for the Williams, MP4exp-6(iii), MP4exp-6(iv), Tsuzuki, OLPS all-atom and MM3 potentials. These correlate fairly well (cc 0.932) with the corresponding parameters of Table 5. The obtained curve leads to the value of 1.158 as the correspondent to the experimental ∆Hvap of 1.96 kcal / mol. The work of Nagy et al.23 presents results of ∆Hvap values for the MM2, MM3, Boyd, MUB-2, Amber, OLPS all-atom, MM2en, MM3en, MM2mc, MM3mc and MM3envir. There is an excellent linear correlation coefficient (c.c. 0.9965) between these ∆Hvap values and the corresponding parameters of Table 5 as shown in Figure 2. The correspondent value of the experimental point is 1.156. Gay et. al.19 carried out Monte Carlo simulations for the E2-Gay, E4-Gay and Williams potentials. These 3 points also correlate with the results of table 5 and leads to the value of 1.152 as the correspondent of the experimental ∆Hvap.


20

Page 21 of 31

MM2en

MM3en

MM2mc

MM3mc

Rowley-A

Rowley-BB

RMK

OLPS allatom

Ouyang

MP4exp6(iii)

MP4exp-6(iv)

Chao

Tsuzuki

MM2

MM3

E2-Gay

E4-Gay

Boyd

Amber

MUB-2

Williams

Sheikh

MG

TraPPE-EH

MM3envir

Table 4: Ratios of the Aggregate Cohesive Energy by the SF Energya

CLC

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48


CLC

0.894

0.882

0.888

0.884

0.885

0.890

0.874

0.833

0.871

0.877

0.874

0.885

0.894

0.887

0.895

0.895

0.896

0.637

0.882

0.887

0.883

0.888

0.894

0.898

0.829

0.873

MM2en

0.947

0.941

0.946

0.946

0.947

0.950

0.944

0.903

0.933

0.932

0.943

0.948

0.947

0.949

0.950

0.952

0.950

0.790

0.944

0.945

0.947

0.950

0.949

0.951

0.932

0.944

MM3en

0.952

0.947

0.951

0.952

0.953

0.956

0.951

0.910

0.940

0.938

0.949

0.954

0.952

0.955

0.956

0.958

0.955

0.804

0.950

0.951

0.953

0.956

0.954

0.956

0.939

0.951

MM2mcb

1.157

1.151

1.155

1.159

1.161

1.161

1.162

1.105

1.151

1.143

1.157

1.160

1.157

1.162

1.163

1.164

1.158

1.001

1.155

1.159

1.162

1.164

1.159

1.162

1.153

1.162

MM3mc

1.153

1.147

1.151

1.154

1.157

1.156

1.157

1.105

1.147

1.139

1.153

1.156

1.153

1.157

1.158

1.159

1.153

1.012

1.151

1.154

1.157

1.159

1.155

1.157

1.150

1.158

Rowley-A

0.823

0.819

0.824

0.824

0.826

0.829

0.822

0.778

0.808

0.806

0.821

0.827

0.823

0.828

0.828

0.831

0.828

0.659

0.822

0.823

0.826

0.828

0.826

0.829

0.809

0.822

Rowley-B

1.138

1.131

1.135

1.142

1.146

1.144

1.150

1.072

1.136

1.122

1.141

1.144

1.138

1.146

1.146

1.147

1.138

0.957

1.137

1.143

1.147

1.149

1.141

1.144

1.142

1.149

RMK

1.206

1.205

1.205

1.198

1.196

1.201

1.188

1.220

1.194

1.213

1.197

1.199

1.206

1.198

1.202

1.201

1.210

1.193

1.202

1.200

1.195

1.196

1.208

1.206

1.178

1.189

OLPS

1.193

1.187

1.189

1.192

1.193

1.192

1.193

1.153

1.191

1.187

1.190

1.192

1.193

1.192

1.195

1.195

1.192

1.071

1.189

1.193

1.193

1.194

1.194

1.195

1.182

1.192

Ouyang

0.970

0.965

0.966

0.963

0.962

0.965

0.955

0.953

0.960

0.974

0.960

0.963

0.970

0.963

0.968

0.967

0.973

0.872

0.964

0.966

0.961

0.963

0.973

0.972

0.939

0.956

POT.

MP4exp-6(iii)

1.191

1.185

1.191

1.194

1.197

1.198

1.197

1.135

1.180

1.172

1.193

1.197

1.191

1.199

1.198

1.200

1.193

1.019

1.190

1.192

1.198

1.200

1.193

1.196

1.191

1.198

MP4exp-6(iv))

1.156

1.151

1.156

1.158

1.160

1.162

1.158

1.105

1.143

1.138

1.156

1.161

1.156

1.162

1.162

1.164

1.159

0.987

1.154

1.156

1.160

1.163

1.158

1.161

1.150

1.159

Chao

1.210

1.201

1.204

1.202

1.202

1.205

1.196

1.172

1.196

1.201

1.196

1.202

1.210

1.204

1.209

1.209

1.210

1.051

1.201

1.205

1.201

1.204

1.209

1.211

1.165

1.195

Tsuzuki

1.167

1.162

1.168

1.171

1.174

1.175

1.174

1.109

1.155

1.146

1.169

1.175

1.167

1.176

1.174

1.178

1.170

0.986

1.167

1.169

1.175

1.177

1.170

1.173

1.167

1.175

MM2

1.772

1.763

1.767

1.770

1.771

1.772

1.769

1.719

1.763

1.760

1.766

1.770

1.772

1.772

1.776

1.776

1.772

1.597

1.766

1.771

1.771

1.774

1.774

1.776

1.754

1.769

MM3

1.008

1.001

1.007

1.009

1.011

1.013

1.010

0.952

0.996

0.990

1.006

1.011

1.008

1.013

1.014

1.016

1.010

0.827

1.005

1.008

1.012

1.014

1.010

1.013

0.998

1.010

E2-Gay

0.954

0.951

0.953

0.950

0.949

0.953

0.943

0.940

0.940

0.950

0.947

0.951

0.954

0.951

0.954

0.954

0.958

0.867

0.951

0.950

0.948

0.951

0.956

0.957

0.929

0.943

E4-Gay

1.149

1.148

1.148

1.147

1.146

1.146

1.145

1.155

1.148

1.152

1.146

1.146

1.149

1.146

1.147

1.146

1.149

1.162

1.147

1.147

1.145

1.146

1.149

1.148

1.143

1.145

Boyd

1.130

1.123

1.128

1.129

1.131

1.134

1.128

1.080

1.116

1.114

1.126

1.131

1.130

1.133

1.134

1.136

1.133

0.957

1.127

1.129

1.131

1.134

1.132

1.135

1.113

1.128

Amber

1.129

1.123

1.126

1.129

1.131

1.131

1.130

1.082

1.125

1.120

1.127

1.130

1.129

1.131

1.133

1.134

1.130

0.978

1.126

1.130

1.131

1.133

1.131

1.133

1.117

1.130

MUB-2

1.199

1.193

1.199

1.203

1.206

1.206

1.207

1.140

1.190

1.180

1.202

1.206

1.199

1.207

1.206

1.209

1.201

1.024

1.198

1.201

1.207

1.209

1.202

1.205

1.202

1.208

Williams

1.133

1.127

1.131

1.134

1.136

1.137

1.136

1.087

1.125

1.120

1.132

1.136

1.133

1.137

1.138

1.139

1.134

0.990

1.131

1.134

1.137

1.139

1.135

1.137

1.128

1.136

Sheikh

1.189

1.178

1.183

1.180

1.180

1.186

1.170

1.132

1.168

1.181

1.174

1.182

1.189

1.183

1.190

1.190

1.194

0.941

1.180

1.183

1.179

1.184

1.193

1.194

1.135

1.170

MG

1.208

1.198

1.204

1.205

1.208

1.211

1.203

1.133

1.192

1.190

1.199

1.207

1.208

1.209

1.213

1.215

1.211

0.942

1.202

1.207

1.208

1.211

1.211

1.215

1.173

1.202

TraPPE-EH

1.103

1.100

1.103

1.107

1.109

1.108

1.111

1.067

1.101

1.092

1.108

1.110

1.103

1.110

1.109

1.111

1.104

1.010

1.104

1.105

1.111

1.112

1.105

1.107

1.115

1.113

MM3envir

1.221

1.215

1.219

1.223

1.226

1.226

1.227

1.168

1.215

1.206

1.221

1.225

1.221

1.226

1.227

1.228

1.222

1.067

1.219

1.223

1.227

1.228

1.223

1.226

1.219

1.227

a: Energy of the aggregate calculated by SF potential is -680.25 kcal/mol b: in bold the ten best performing potentials 21 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 31

Figure 2: Linear correlation between the Aggregate energy ratio and the per cent deviation of ∆Hvap obtained by MC simulation of 11 different potentials by Nagy et al23.

The value of 1.156 was taken as the reference by which the potentials studied are measured, i.e., the potentials that present the closest value to 1.156 are considered to provide the best cohesive energies. These correlations show that the aggregate energies can be taken as rough approximate surrogate values to the latent heat of evaporation of methane (∆Hvap). It is helpful in a study like this one in which a large number of potentials are considered and for many of them no ∆Hvap has been reported in the original literature. The ten closest values to 1.156 at Table 5 belong to the series: MP4exp-6(iv), MM2mc, MM3mc, E4-Gay, Tsuzuki, Rowley-B, Williams, Sheikh, Amber and Boyd.


22

Page 23 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


A Measurement of the Anisotropy of the Potentials. The proposed scale of anisotropy (Eopt – Erand)/Eopt where Eopt is the energy of the optimized and Erand is the energy of the randomized aggregate is presented in Table 5. The smaller these numbers are the less important are the relative positions of the hydrogen atoms for each potential. At one extreme lie the united-atom models with a ratio of zero and at the other end lies potential CLC with a value of 2.27. The latter is a four site model where the carbon atom is not an interacting site. One may speak now of a “united-atom character” of a potential to designate how well it resembles the united atom potentials or alternatively it can be understood as a measure of the anisotropy of the potentials.

Table 5: Anisotropy, Ratio of the Aggregate Energy, RMS of the Dimers Energy for the Data of Chao-Li, Tsuzuki and Hayes POTENTIAL ANISOTROPYa Aggregate RMSb dimer RMSb RMSb energy Chao-Li dimer dimer Tsuzuki Hayes CLC MM2enc MM3en MM2mc MM3mc Rowley-A Rowley-B RMK OLPS all-atom Ouyang MP4exp-6(iii) MP4exp-6(iv) Chao Tsuzuki MM2 MM3 E2-Gay E4-Gay SF Boyd Amber MUB-2

2.272 0.370 0.341 0.277 0.239 0.396 0.353 0.121 0.471 0.584 0.246 0.279 1.130 0.260 0.273 0.366 0.309 0.442 0.000 0.349 0.558 0.232

0.884 ± 0.014 0.945 ± 0.006 0.951 ± 0.005 1.158 ± 0.005 1.154 ± 0.005 0.823 ± 0.007 1.142 ± 0.006 1.200 ± 0.008 1.192 ± 0.003 0.964 ± 0.007 1.193 ±0.007 1.157 ± 0.006 1.202 ± 0.009 1.170 ± 0.007 1.769 ± 0.005 1.008 ± 0.006 0.950 ± 0.006 1.147 ± 0.002 1.000 1.128 ± 0.006 1.129 ± 0.004 1.202 ± 0.007

0.225 0.064 0.056 0.058 0.051 0.077 0.090 0.082 0.087 0.079 0.051 0.044 0.109 0.056 0.150 0.064 0.105 0.164 0.054 0.096 0.062


2.455 0.150 0.116 0.139 0.111 0.150 0.292 0.429 0.555 0.565 0.106 0.090 1.578 0.100 0.338 0.145 1.977 2.821

7.656 0.193 0.157 0.248 0.172 0.160 0.681 0.370 1.956 1.885 0.154 0.138 5.148 0.150 0.512 0.222 1.892 4.715

0.186 0.639 0.181

0.278 2.143 0.228

23


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Williams Sheikh HFD MG TraPPE-EH MM3envir

0.247 1.150 0.000 1.118 0.124 0.241

1.134 ± 0.005 1.181 ± 0.012 1.033 1.205 ± 0.009 1.107 ± 0.005 1.223 ± 0.005

0.050 0.210 0.191 0.084 0.012

Page 24 of 31

0.137 1.643

0.193 4.896

1.574 0.935 0.132

4.820 1.015 0.194

a: measured as (Eopt – Erand)/Eopt ,b: in kcal/mol, c: in bold the ten best performing potentials as regards dimer energies. DISCUSSION Each one of the three studied properties yielded a set of the ten best performing potentials but only three potentials are present simultaneously at the three sets: MM3mc, MM2mc and Tsuzuki. The set of the top ten multi-site potentials based on the second virial coefficients have two points in common: 1 - They all have larger aggregate cohesive energies compared to the reference value of 1.156 which indicates that these potentials overestimate the heat of evaporation and 2 - they have smaller anisotropies (see Table 5). The latter finding is consistent with the observation that the united-atom models have generally better performances at the second virial coefficients. The potential TraPPE-EH has a lower aggregate cohesive energy (1.107) but in this case it is likely to show a slightly higher ∆Hvap than the experimental value because it shows a heavy united-atom character (anisotropy equals 0.124, Table 5). Once again, the notable exception is E4-Gay, which is the second best performer at the virials and its value at Table 5 is 0.442 showing a large anisotropy. It also does very well with the aggregate cohesive energies which corroborates the good ∆Hvap obtained by Gay et.al.19 but does very poorly at the reproduction of dimers ab initio interaction energies. In fact, potential E4-Gay is so unlike every other potential including ab initio data that it should be seen with caution.


24

Page 25 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The energies obtained with the aggregates show that the studied potentials are very similar in recognizing the most optimized structure. This is consistent with their agreement observed in the order of increasing r*, the location of the energy minimum of the orientations of the dimers. Curiously, even effective potentials such as OLPS allatom, Amber, Ouyang, Chao, Boyd, Williams agree with the ab initio data. The top ten potentials based on the aggregate energies differ significantly from those obtained with the second virial coefficients. This apparent contradiction may be understood on the basis of the non-additivity of the potentials. Following the view of Hauschild and Prausnitz33 and of Abbaspour24, who have explicitly taken three body effects in consideration along with two body potentials, the many body potentials are positive and diminishes the total cohesiveness of the aggregates. Thus, effective potentials, constructed to model the bulk properties of liquids will do so nicely but in comparison to true pair potentials will be less attractive. Therefore, it is not reasonable to expect that a potential will yield simultaneously accurate second virial coefficients and good latent heat of evaporation unless some form of many body correction is employed. This interpretation is in perfect agreement with the fact that so many potentials obtained in a variety of ways cannot reproduce these two quantities simultaneously. Except for RMK and E4-Gay already discussed, and the known deviation of MM2, all potentials underestimate the second virial coefficients. Szczesniak et al.11 conducted ab initio calculation of the dispersion energy in trimers of methane and concluded that the three body potential is very small. But these authors have indeed obtained shallow potential wells even for the dimers compared to later ab

initio computations. Thus it is fair to expect that they have also greatly underestimated the extent of the three body potential.


25


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 31

The potentials obtained through ab initio calculations of the dispersion energy of dimers: MM3en and MM2en of Nagy et. al.18 , MP4exp-6(iii) and MP4exp-6(iv) of Hayes et al10, Tsuzuki, CLC, Chao, Rowley-A and Rowley-B should in principle yield accurate second virial coefficients rather than good latent heat of evaporation because energy of dimers are related to very low density where many body effects are absent. The potentials Tsuzuki and MP4exp-6(iii), obtained solely by ab initio methods, have yielded very good results at the second virial coefficients and are ranked 5th and 6th among the top ten multi-site potentials but still show B(T) values approximately 10-15 % lower than the experimental ones. This indicates that ab initio methods have not yet attained chemical accuracy as far as potential well depths due to dispersion energies are concerned. Tsuzuki et al.34. in 2006 revisited the calculation of the potential well depth for the J orientation of the methane dimer and found it to be even deeper, at -0.48 kcal/mol, about 9 % deeper than previous values. The picture that emerges from the above discussion is that the best pair potential is the MM3envir of Nagy et al.23 which possess anisotropy of 0.241 but can discern very well the orientations of the dimers, followed closely by the united atom potential of Abbaspour24 (HFD). MUB-2 is also a very interesting potential yielding good second virial coefficients and reasonable dimers energy. The RMK potential is the best at reproducing second virial coefficients but has a very small anisotropy (0.121) and likewise E4-Gay disagrees very much with ab initio data. For the condensed phase, the best potentials would be those that appear simultaneously at the top ten of the aggregate energies and dimers energy: MP4exp-6(iv), MM2mc, MM3mc, Tsuzuki and Williams. Interestingly, in this lot there are two potentials obtained purely from ab initio data: Tsuzuki and MP4exp-6(iv), two obtained


26

Page 27 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


simultaneously from bulk properties and ab initio data: MM2mc and MM3mc and one obtained from crystal lattice energies- Williams. Their anisotropies lie in the range 0.239 to 0.279. The strategy of Nagy et. al.23 of parameterizing the MM2mc and MM3mc potentials with more than one property seems to have been effective. These potentials also happen to be in the top ten of second virial coefficients. CONCLUSION

A molecular dynamics simulation of liquid methane at its normal boiling point using the Saager Fischer (SF) potential yielded the carbon infrastructure needed to build aggregates of 512 molecules. As their orientation were optimized by each studied potential an approximate surrogate parameter for the latent heat of evaporation emerged. Based on this parameter, the aggregate cohesive energy, a ranking of the studied multisite potentials was established. Though only approximate, this surrogate parameter is useful to select those potentials that yield good ∆Hvap values without having to resort to MD simulations to all studied potentials. It was found that the best potentials to use in condensed phases are: MP4exp-6(iv) of Hayes et al.10, MM2mc and MM3mc of Nagy

et. al.23, the potential of Tsuzuki21 and the potential of Williams1. The optimized aggregates were common to all potentials, excepting E4-Gay and RMK ones. This finding is consistent with the behavior of most potentials in following the same order of increasing r*, the location of the energy minimum for each orientation of the dimers. By analyzing the difference in energy between the fully optimized and the randomly oriented aggregates it was possible to set up a quantitative measurement of the anisotropy of these potentials.


27


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 31

The potentials with the smaller anisotropies showed to be better in reproducing second virial coefficients. The MM3envir potential of Nagy et. al23. was found to be the closest one to a true pair potential. Potentials RMK and E4-Gay are respectively the first and second in the ranking of the virial coefficients but they don’t reproduce the ab initio behavior. The potentials that best reproduce second virial coefficients tend to overestimate the latent heat of evaporation despite not being attractive enough to totally satisfy the virials. This can be understood on the basis of the non-additivity of the potentials as pointed out by Hauschild and Prausnitz33. The non-additivity will make a pair potential increasingly less attractive as density increases due to many body interactions. As even the best potentials derived from ab initio data missed the experimental values of second virial coefficients by 10-15% it is likely that these methods have not yet attained chemical accuracy as far as calculating dispersion energy is concerned.

ACKNOWLEDGEMENTS The author is grateful to Dr. Alejandro Lopes Castillo of the Federal University of São Carlos – Brazil and to Dr. Luzia Peres Novaki of the University of São Paulo – Brazil for commenting on the manuscript prior to submission.

REFERENCES 1- Williams, D.E. Nonbonded Potential Parameters Derived from Crystalline Hydrocarbons. J. Chem. Phys., 1967, 47, 4680-4684. 2- Aimoli, C.G.; Maginn, E.J.; Abreu, C.R.A. Force Field Comparison and Thermodynamic Property Calculation of Supercritical CO2 and CH4 Using Molecular Dynamics Simulations. Fluid Phase Equilibria. 2014, 368, 80-90. 3- Ungerer, P.; Lachet, V.; Tavitian, B. Applications of Molecular Simulation in Oil and Gas Production and Processing, Oil & Gas Science and Technology-Rev.IFP, 2006, 61, 387-403


28

Page 29 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


4- Younglove, B.A.; Ely, J.F. Thermophysical Properties of Fluids. II. Methane, Ethane, Propane, Isobutane, and Normal Butane. J. Phys. Chem. Ref. Data, 1987, 16, 577-798. 5- Stassen, H. On the Pair Potential in Dense Fluid Methane. J. Mol. Struct. (Theochem), 1999, 464, 107-119. 6- Skarmoutsos, I.; Kampanakis, L.I.; Samios, J. Investigation of the Vapor-Liquid Equilibrium and Supercritical Phase of Pure Methane via Computer Simulations. J. Mol. Liquids, 2005, 117, 33-41. 7- Saager, B.; Fischer, J. Predictive Power of Effective Intermolecular Pair Potentials: MD Simulation Results for Methane up to 1000 MPa. Fluid Phase Equilibria, 1990, 57, 35-46. 8- Jorgensen, W.L.; Madura, J.D.; Swenson, C.J. Optimized Intermolecular Potential Functions for Liquid Hydrocarbons. J. Am. Chem. Soc, 1984, 106, 6638-6646. 9- Chao, S.W.; Li, A.H.T.; Chao, S.D. Molecular Dynamics Simulations of Fluid Methane Properties Using Ab Initio Intermolecular Interaction Potentials. J. Comput. Chem., 2008, 30, 1840-1849. 10- Hayes, J.M.; Greer, J.C.; Morton-Blake, D.A. A Force Field Description of ShortRange Repulsions for High Alkane Molecular Dynamics Simulations. J. Comput. Chem., 2004, 25, 1953-1966. 11- Szczesniak, M.M.; Chafasinski, G.; Cylbulski, S.M.; Scheiner, S. Intermolecular Potential of the Methane Dimer and Trimer., J. Chem. Phys. 1990, 93, 42434253. 12- Chen, B.; Siepmann, J.I. Transferable Potentials for Phase Equilibria. 3. ExplicitHydrogen Description of Normal Alkanes. J. Phys. Chem. B. 1999, 103, 53705379. 13- Kaminski, G.; Duffy, E.M.; Matsui, T.; Jorgensen, W.L. Free Energies of Hydration and Pure Liquid Properties of Hydrocarbons from the OLPS All-Atom Model. J. Phys. Chem., 1994, 98, 13077-13082. 14- Allen, M.P.; Tildesley, D.J. Computer Simulation of Liquids, Clarendon Press – Oxford, 1987 15- Lii, J.-H.; Allinger, N.L. Molecular Mechanics. The MM3 Force Field for Hydrocarbons. 3. The Van der Waals’ Potentials and Crystal Data for Aliphatic and Aromatic Hydrocarbons. J. Am. Chem. Soc. 1989,111, 8576-8582. 16- Righini, R.; Maki, K.; Klein, M.L. An Intermolecular Potential for Methane, Chem. Phys. Lett., 1981, 80, 301-305. 17- Allinger, N.L. Conformational Analysis. 130. MM2. A Hydrocarbon Force Field Utilizing V1 and V2 Torsional Terms. J. Am. Chem. Soc, 1977, 99, 8127-8134. 18- Fitzwater, S.; Bartell, L.S. Representations of Molecular Force Fields. 2. A Modified Urey-Bradley Field and an Examination of Allinger’s gauche Hydrogen Hypothesis. J. Am. Chem. Soc, 1976, 98, 5107-5115.


29


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 30 of 31

19- Gay, D.H.; Dai, H.; Beck, D.R. Obtaining Accurate Pressure Second Virial Coefficients for Methane from Ab Initio Pair Potential. J. Chem. Phys., 1991, 95, 9106-9114. 20- Jiang, C.; Ouyang, J.; Zhuang, X.; Wang, L.; Li, W. An Efficient Fully Atomistic Potential Model for Dense Fluid Methane. J. Mol. Struct., 2016, 1117, 192-200. 21- Tsuzuki, S.; Uchimaru, T.; Tanabe, K. A New Ab Initio Based Model Potential for Methane., Chem. Phys. Lett., 1998, 287, 327-332. 22- Rowley, R.L.; Pakkanen, T.A. Determination of a Methane Intermolecular Potential Model for Use in Molecular Simulations from Ab Initio Calculations. J. Chem. Phys. 1999, 110, 3368-3377. 23- Nagy, J.; Weaver, D.F.; Smith, V. H. Jr. A Comprehensive Study of Alkane Nonbonded Empirical Force Fields. Suggestions for Improved Parameter Sets. J. Phys. Chem., 1995, 99, 8058-8065. 24- Abbaspour, M. Computation of Some Thermodynamics, Transport, Structural Properties, and New Equation of State for Fluid Methane Using Two-Body and Three-Body Intermolecular Potentials from Molecular Dynamics Simulation. J. Mol. Liquids, 2011, 161, 30-35. 25- Li, A.H.T.; Chao, S.D. A Refined Intermolecular Interaction Potential for Methane: Spectral Analysis and Molecular Dynamics Simulations. J. Chin. Chem. Soc., 2016, 63, 282-289. 26- Chang, S.-J.; McNally, D.; Shary-Tehrani, S.; Hickey, M.J.; Boyd, R.H. Heats of Combustion and Strain Energy of Bicyclo[n.m.O]alkanes. J. Am. Chem. Soc, 1970, 92, 3109-3118. 27- Sun, Y.; Spellmayer, D.; Pearlman, D.A.; Kollman, P. Simulation of the Solvation Free Energies for Methane, Ethane, and Propane and Corresponding Amino Acid Dipeptiedes: A Critical Test of the “Bond-PMF” Correction, a New Set of Hydrocarbon Parameters, and the Gas Phase-Water Hydrophobicity Scale. J. Am. Chem. Soc, 1992,114, 6798-6801 28- El-Sheikh, S.M.; Barakat, K.; Salem, N.M. Phase Transitions of Methane Using Molecular Dynamics Simulations. J. Chem. Phys. 2006, 124, 124517 29- Murad, S.; Gubbins, K.E. Molecular Dynamics Simulation of Methane Using a Singularity-Free Algorithm. ACS Symp. Ser., 1978, 86, 62-71. 30- McQuarrie, D.A. Statistical Mechanics, Harper and Row, New York, 1976 31- Dymond, J.H.;Smith, E.B. The Virial Coefficients of Pure Gases and Mixtures, Clarendon, Oxford, 1980. 32- Rapaport, D.C. The Art of Molecular Dynamics Simulation, Cambridge University Press, Cambridge, 1995 33- Hauschild, T.; Prausnitz, J.M. Monte Carlo Calculations for Methane and Argon Over a Wide Range of Density and Temperature, Including the Two-Phase Vapor-Liquid Region., Mol. Simul., 1993, 11, 177-185.


30

Page 31 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


34- Tsuzuki, S.; Honda, K.; Uchimaru, T.; Mikami, M. Estimated MP2 and CCSD(T) Interaction Energies of n-Alkane Dimers at the Basis Set Limit: Comparison of the Methods of Helgaker et al. and Feller., J. Chem. Phys., 2006, 124, 114304.


31

Intermolecular Potentials of Methane Assessed by Second Virial

Recommend Documents