Calculating the Fugacity of Pure, Low Volatile Liquids via Molecular

Article pubs.acs.org/IECR

Calculating the Fugacity of Pure, Low Volatile Liquids via Molecular Simulation with Application to Acetanilide, Acetaminophen, and Phenacetin Georgia B. Fuerst,† Ryan T. Ley,† and Andrew S. Paluch* Department of Chemical, Paper and Biomedical Engineering, Miami University, Oxford, Ohio 45056, United States

Downloaded by SUNY UPSTATE MEDICAL UNIV on September 5, 2015 | http://pubs.acs.org Publication Date (Web): September 1, 2015 | doi: 10.1021/acs.iecr.5b01827

S Supporting Information *

ABSTRACT: Conventional molecular simulation free energy calculations and standard thermodynamic relations are applied to compute the pure liquid fugacity of low volatile liquids and compared to reference Monte Carlo simulations. The method involves the calculation of the residual chemical potential and the molar volume of the liquid at the conditions of interest. For substances that are solid at the conditions of interest, simulations may be performed at elevated temperatures and extrapolated to subcooled conditions. It is shown that direct calculations at subcooled conditions provide erroneous values of the fugacity. Knowledge of the pure liquid fugacity is essential to compute activity coefficients defined with respect to a Lewis−Randall standard state for thermodynamic property modeling. The validity of the method is verified by comparing to reference simulation results for n-octane, 3,4-dimethylhexane, cyclohexane, methanol, 1-propanol, and 2-propanol. The method is then applied to acetanilide, acetaminophen, and phenacetin, all of which are solid at ambient conditions. The results for the six reference compounds are in good agreement with standard Monte Carlo simulations, suggesting that free energy based calculations of the fugacity may be used to accurately normalize activity coefficients of low volatile liquids computed via molecular simulation.

1. INTRODUCTION A quantitative description of the phase behavior of molecular systems is crucial for the development of efficient separation processes. Historically, the necessary thermodynamic data are obtained from empirical or semiempirical correlations of experimental data.1,2 However, it is challenging to apply these conventional approaches to novel compounds encountered during design processes. Early in the design process, the necessary data likely do not exist and may be challenging to obtain. It may be the case that the proposed compound may not yet have been synthesized. It follows that accurate predictive methods would be of great utility.3−7 The present study will focus on the calculation of the pure liquid fugacity of low volatile liquids. The ability to accurately determine the pure liquid fugacity is primarily desirable because it is the normalization constant necessary to calculate the activity coefficient of a nonelectrolyte solid solute in solution defined with respect to a Lewis−Randall (or Raoult’s law) standard state.8,9 For liquids at low pressure, the pure liquid fugacity is typically taken to be equal to the vapor pressure of the pure liquid. Knowledge of the vapor pressure of (low volatile) pesticides is required by the U.S. Environmental Protection Agency.10 For nonelectrolyte compounds that are solid at ambient conditions, the fugacity is needed for a subcooled liquid state. For this case, the vapor pressure is measured at elevated temperatures and extrapolated down to the temperature of interest. However, accurately measuring the vapor pressure of low volatile compounds experimentally is with great challenge, requiring specialized techniques.10−15 This challenge is further complicated for species of pharmaceutical interest where the range of stability of the compound in vapor− liquid equilibrium is limited. © XXXX American Chemical Society

Here we will demonstrate how conventional molecular simulation free energy calculations may be used to accomplish this task. Motivated by the desire to accurately compute the pure liquid fugacity of nonelectrolyte compounds that are solid at ambient conditions, the method used in the present study is particularly well suited for compounds below their normal boiling point but above their normal melting point. When used in combination with free energy calculations for the compound in a solvent of interest, the activity coefficient of the compound may be computed. This will allow for the use of established solution theories and excess Gibbs free energy models in combination with molecular simulation to quantitatively predict solid−liquid phase behavior, in addition to allowing new solution theories and excess Gibbs free energy models to be developed. This will be of great value, for instance, for solvent selection during the design of extraction and crystallization processes involving novel active pharmaceutical ingredients,5 the extraction of pollutants from waste- and groundwater,16,17 and the desulfurization of diesel fuel.18 Free energy based calculations of the fugacity are appealing, as they use the same free energy based molecular simulation methodology used to compute solution phase free energies, which is implemented in many molecular simulation software (both molecular dynamics and Monte Carlo), thereby facilitating the adoption of the method by novice practitioners. Additionally, use of the method to estimate the vapor pressure of low volatile liquids is appealing, for which experimentation Received: May 17, 2015 Revised: August 14, 2015 Accepted: August 26, 2015

A

DOI: 10.1021/acs.iecr.5b01827 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article


Industrial & Engineering Chemistry Research and conventional molecular simulation methods are challenging. The method is first applied to n-octane, cyclohexane, 3,4dimethylhexane, methanol, ethanol, and 2-propanol for which ample reference simulation data exist, in addition to reference Monte Carlo simulation results generated in the present study. The method is then applied to acetanilide, acetaminophen, and phenacetin, all of which are solid at ambient conditions. Calculations are performed between the experimental normal boiling and melting point, which may be extrapolated to subcooled conditions. It is additionally demonstrated with acetaminophen that direct calculations below the normal melting point yield erroneous results. Although not utilized in the present study, in the Supporting Information accompanying this manuscript we additionally show how results generated from fugacity calculations may be used to compute the Hildebrand solubility parameter. This is useful for efficiently predicting solubility trends with regular solution theory.8,9

which is in agreement with the expression derived by Lin et al.22 Note that the units of f L are the same as RT/v (i.e., units of pressure); throughout the present study, f L will have units of kPa unless noted otherwise. For large N, eq 5 is rigorously correct and is the basis of the present study. If the vapor phase in equilibrium with the liquid phase at T is assumed to be an ideal gas (such that the fugacity coefficient is unity) and the Poynting correction is negligible, then f L ≈ psat, where psat is the saturation pressure. For this case, eq 5 reduces to the expression derived in the work of Winget et al.23 who derived an expression to calculate the vapor pressure of liquids at 298 K using solvation free energy calculations assuming “ideal behavior.” In the Supporting Information accompanying this manuscript, we derive the expression for psat from Winget et al. in an analogous fashion as the NpT plus test particle method of Fischer and co-workers24−26 and show that ideal behavior corresponds to assuming that the liquid phase chemical potential is independent of pressure and that the vapor phase is an ideal gas. With these two assumptions, the same expression was derived by Horn et al.27 who cleverly used it to predict the vapor pressure and normal boiling point of their novel TIP4P-Ew water model using molecular dynamics simulations in an NpT ensemble. Additionally, the expression for psat from Winget et al. was recently used by Ahmed and Sandler28 to predict the subcooled liquid vapor pressure of nitroaromatic energetic compounds at ambient conditions to ultimately predict their equilibrium solubility. The systems of interest for drug design are solid at ambient conditions. Therefore, eq 5 would require performing simulations of a subcooled liquid. While this could be done as in ref 28, extreme care must be taken when modeling subcooled liquids. Such systems are quasi-nonergodic, requiring long simulation times and/or sophisticated enhanced sampling techniques to obtain accurate thermodynamic properties.29−32 The risk of performing free energy calculations for subcooled liquids will be further emphasized in the “Results and Discussion” section. For subcooled liquids, the results of the present study suggest it is better to calculate f L for a range of temperatures above the normal melting point then extrapolate the results to ambient conditions. 2.2. Grand-Canonical Transition-Matrix Monte Carlo. The fugacity at saturation may be rigorously computed using grand-canonical transition-matrix Monte Carlo (GC-TMMC) simulation33−35 with histogram reweighting.26,32,36,37 The methodology is well established and will only be summarized here. The interested reader can find an expanded discussion in the Supporting Information accompanying this manuscript and detailed discussions in refs 32−35. In a GC-TMMC simulation, one specifies the system temperature, volume, and chemical potential and computes the molecule number probability distribution. After the simulation is complete, the probability distribution may be obtained at another chemical potential with use of histogram reweighting. The chemical potential at saturation is found by reweighting the probability distribution to determine the value of the chemical potential that yields a bimodal probability distribution. Once the probability distribution at saturation is known, the saturation pressure and molar volume of the liquid and vapor phase may readily be computed.33,37 The fugacity at saturation may then be computed as38,39

2. METHODOLOGY 2.1. Free Energy Calculations. Performing conventional free energy calculations in the isothermal−isobaric (NpT) ensemble, one may readily calculate the residual chemical potential (or the residual molar Gibbs free energy) of a pure liquid as19−21 βμres (T ,p) = βμL (T ,p) − βμig (T ,p′)

(1)

where T is the temperature, β−1 = kBT where kB is Boltzmann’s constant, p is the pressure, μres and μL are the residual chemical potential and the liquid phase chemical potential of the pure liquid at T and p, respectively, and μig is the chemical potential of the species in a noninteracting, ideal gas state at the same T but at pressure p′, where NkBT ⟨V ⟩T , p , N − 1

p′ =

(2)

where N is the number of molecules in the system and ⟨V⟩T,p,N−1 is the ensemble average volume of the system in the absence of the molecule that is being coupled/decoupled to the system.19−21 By use of fundamental thermodynamic relations, μres may be related to the liquid phase fugacity, f L, as8 f L (T ,p) = ln f L (T ,p) − ln p′ f ig (T ,p′) NkBT = ln f L (T ,p) + ln ⟨V ⟩T , p , N − 1

βμres (T ,p) = ln

(3)

ig

where f is the fugacity of the species in the noninteracting, ideal gas state (i.e., the same conditions as μig in eq 1). Next, assume that N ≫ 1 such that ⟨V ⟩T , p , N − 1 NkBT

≈

⟨V ⟩T , p , N NkBT

=

v(T ,p) RT

(4)

where v is the liquid molar volume and R = NavokB is the universal gas constant where Navo is Avogadro’s number. The final expression for the pure liquid fugacity is given as ln f L (T ,p) = βμres (T ,p) + ln

RT v(T ,p)

⎛k T ⎞ ln f sat = βμsat + ln⎜ B3 ⎟ + ln Zintra ⎝ Λ ⎠

(5) B

(6)


Article

Industrial & Engineering Chemistry Research


where Λ = h/(2πmkBT)1/2 is the thermal de Broglie wavelength, h is Planck’s constant, m is the mass of a molecule, and Zintra is the (intramolecular) ideal gas configurational integral of a single molecule. Zintra is only dependent on T and may readily be evaluated via Monte Carlo integration by performing a simulation with a single molecule. For molecular species, Zintra is calculated as the ensemble average of the Rosenbluth weight computed for growing a single molecule in vacuum (i.e., in an isolated ideal gas state).40 For species that lack any intramolecular degrees of freedom (such as a LennardJones fluid), Zintra = 1. For liquids at temperatures below the normal boiling point (which are the motivation for the present study), the fugacity coefficient at saturation and the Poynting correction are expected to be near unity such that the vapor pressure, fugacity at saturation, and the liquid fugacity may be assumed numerically equivalent.8

pseudoatom with TraPPE-UA LJ parameters.41,51 The LJ parameters for the hydroxyl group of acetaminophen were taken from phenol in TraPPE-EH,47 and the ethoxy group LJ parameters for phenacetin were taken from the TraPPE-UA ether model.41,50,51 Partial atomic charges were obtained for the TraPPE-EH based models using a protocol similar to the original TraPPEEH parametrization.45−47 Slight differences emerged due to the capabilities of the ab initio quantum chemistry software package available to us at the time of this study and the use of unitedatom alkyl sites. Charges were parametrized following a threestep procedure. First, the gas-phase structure for each molecule was optimized at the M06-2X/cc-pVTZ level of theory/basis set.52,53 Second, a single point energy calculation was performed on the gas-phase optimized structure at the M062X/6-31G(d) level of theory/basis set52,53 in a self-consistent reaction field (SCRF) using the SMD universal solvation model for 1-octanol.54−56 All of the electronic structure calculations were performed using Gaussian 09.57 Third, partial atomic charges were then obtained from the electrostatic potential (obtained in step 2) using the restrained electrostatic potential (RESP)58,59 method in ANTECHAMBER (part of the AMBER 12 simulation suite).60,61 During the fitting procedure, the partial charges of alkyl hydrogens were fixed to zero; this had the effect of using a single charge site at the center of the alkyl group appropriate for use with a united-atom model. All of the intramolecular parameters for acetanilide, acetaminophen, and phenacetin were taken from the general AMBER force field (GAFF).62,63 Parameters were generated using ANTECHAMBER and converted from AMBER to GROMACS format using ACPYPE.64,65 Since alkyl groups were modeled as a single alkyl pseudoatom, intramolecular parameters involving alkyl hydrogens were removed. Throughout the present study, all bond lengths of all species were held fixed. All of the Gaussian 09 input files and GROMACS force field files for acetanilide, acetaminophen, and phenacetin used in the present study may be found in the Supporting Information of this manuscript. 3.2. Molecular Dynamics. All of the molecular dynamics (MD) simulations were performed with GROMACS 4.6.3.66−68 The MD simulations for each system involved five steps: initial structure generation, minimization, NpT equilibration and production, and free energy calculations. First, initial structures were generated using Packmol.69,70 For n-octane, 3,4dimethylhexane, cyclohexane, methanol, 1-propanol, and 2propanol, the number of molecules was chosen to obtain an equilibrated cubic box with an edge length of approximately 40 Å based on estimates using saturation densities from the original TraPPE-UA publications.41−44 For acetanilide, acetaminophen, and phenacetin, the liquid molar volume at 298 K was estimated using the group contribution method of Hukkerikar et al.7 and used to obtain an estimate of the number of molecules in a cubic box with an edge length of 40 Å. This estimate was used for all simulations at all temperatures. Second, up to 3000 steps of steepest descent minimization were performed on the initial structure to remove any bad contacts that may have been created during packing. For steps 3−5 for dynamics, the equations of motion were integrated with the GROMACS “stochastic dynamics” integrator, corresponding to stochastic or velocity Langevin dynamics integrated with the leapfrog algorithm.71−73 The third step consisted of 100 ps of NpT equilibration with the

3. COMPUTATIONAL DETAILS 3.1. Force Fields. All interactions were modeled using a conventional “class I” potential energy function with all nonbonded intermolecular interactions treated using a combined Lennard-Jones (LJ) and fixed point charge model. The reference compounds n-octane, 3,4-dimethylhexane, cyclohexane, methanol, 1-propanol, and 2-propanol were all modeled using the united atom transferable potentials for phase equilibria (TraPPE-UA) force field.41−44 Molecular models for acetanilide, acetaminophen, and phenacetin (see Figure 1) were constructed based on the

Figure 1. Chemical structures of acetanilide (CAS number 103-84-4), acetaminophen (CAS number 103-90-2), and phenacetin (CAS number 62-44-2).

explicit hydrogen transferable potentials for phase equilibria (TraPPE-EH) force field.45−47 All three molecules have a scaffold consisting of N-phenylacetamide. While acetanilide contains no further functional groups, acetaminophen and phenacetin contain a hydroxyl and ethoxy group attached to the para position, respectively. Beginning with acetanilide, the LJ parameters for the phenyl group were taken directly from TraPPE-EH.45 The LJ parameters for N and H in the secondary amide were taken to be the same as for N and H in a secondary amine in TraPPE-EH.48,49 The carbonyl group LJ parameters were taken from the TraPPE-UA amide model,48,50 and the terminal alkyl (CH3) group was modeled as a single alkyl C


Article



the system (i.e., when λLJ m ≈ 0), intermolecular nonbonded LJ interactions involving the solute were modeled with a modified “soft-core” potential,87−89 and the intermolecular electrostatic interactions are decoupled linearly. Independent MD simulations are performed in each stage m. Periodically throughout these simulations, the change in the Hamiltonian with the current configuration between stage m and every other stage is computed. This information is saved to disk for postsimulation analysis with MBAR86,90−92 to compute μres. In the present study, the solute was taken from noninteracting (m = 0) to fully interacting (m = 9) by first bringing the intermolecular LJ interaction to full strength over 5 stages (1 ≤ m ≤ 5) and then adding in intermolecular electrostatic interactions in the final 4 stages (6 ≤ m ≤ 9) for a total of 10 stages (including the ideal gas reference state m = 0). In stage m = 0, the intermolecular LJ and electrostatic interactions elec between the solute and the system are turned off (λLJ 0 = λ0 = 0). Over the next 5 stages, the intermolecular electrostatic interactions remain off and the intermolecular LJ interactions are strengthened linearly as λLJ m = {0.2, 0.4, 0.6, 0.8, 1.0}. Next, with the LJ interactions fully restored (λmLJ = 1) the intermolecular electrostatic interactions were strengthened following a square root relation as λelec m = {0.50, 0.71, 0.87, 1.00}.81 During this process all bonded and intramolecular LJ and electrostatic interactions are unchanged. The final configuration from the corresponding 15 ns NpT production simulation was taken as the starting configuration in each stage. The NpT ensemble MD simulation in each stage was 10.5 ns in length, with the first 0.5 ns discarded from analysis as equilibration. During these simulations, the differences in the Hamiltonians were computed every 0.20 ps. μres and the associated error was computed using the GROMACS analysis script distributed with the Python implementation of MBAR (PyMBAR).86,90−92 3.3. Monte Carlo. GC-TMMC simulations were additionally performed with n-octane, 3,4-dimethylhexane, and cyclohexane in order to validate the methodology presented in section 2.1. The Monte Carlo (MC) simulations were all performed with the MCCCS Towhee version 7.0.6 simulation package.93,94 Prior to performing the actual GC-TMMC simulations, MC simulations were performed to compute Zintra. For these calculations, the ensemble average Rosenbluth weight was computed using 2 million molecule growths in vacuum. Next, the GC-TMMC simulations were performed. Simulations were performed at a range of temperatures from below the critical point to below the normal boiling point such that the MD and MC simulations have at least one mutual temperature. The lower bound of this range was limited by the feasibility of performing the necessary calculations. The system volume was set to be a cubic box with edge length of 35 Å. In all cases, the LJ interactions were truncated at rcut = 14 Å, with long-range corrections applied to account for the truncation of these interactions.40,79 A hard inner cutoff of rinner cut = 2 Å was also used such that if for two sites (i and j) rij≤ rinner cut , the MC move was immediately rejected or the weight of the CBMC trial site was set to zero. Lorentz−Berthelot combining rules were applied for interactions between unlike LJ sites.79 For the TraPPE-UA models for alkanes, there are no partial charges and all bond lengths are fixed. All GC-TMMC simulations were initiated with an empty box. Trial moves consisted of 15% center-of-mass translations, 15% rotations about the center-of-mass, 15% configurationalbias Monte Carlo (CBMC) regrowths, and 55% CBMC

Berendsen barostat,71,74 followed in step 4 with a 15 ns NpT production simulation with the Parrinello−Rahman barostat.71,75,76 The Parrinello−Rahman barostat was additionally used in step 5 when performing free energy calculations in the NpT ensemble (which will be discussed momentarily). In all three steps, the friction constant for the stochastic (or Langevin) thermostat was 1.0 ps−1, the time constant for the barostat was 5.0 ps, and the pressure was set to 1 bar, and all bond lengths were constrained using P-LINCS.71,77,78 For acetaminophen at 300, 350, and 400 K, steps 4 and 5 were additionally repeated with the friction constant for the stochastic thermostat decreased from 1.0 to 0.125 ps−1. The LJ interactions were switched off smoothly between 13.5 and 14.0 Å, and long-range analytic dispersion corrections were applied to the energy and pressure to account for the truncation of these interactions.40,71,73,79 For interactions between unlike LJ sites Lorentz−Berthelot combining rules were employed.79 Electrostatic interactions were evaluated using the smooth particle-mesh-Ewald (SPME) method with tinfoil boundary conditions71,73,80,81 with real space interactions truncated at 16.5 Å. A SPME B-spline order of 6, a Fourier spacing of 1 Å, and a relative tolerance between long- and short-range energies of 10−6 was used. The molar volume necessary for eq 5 was computed during the 15 ns NpT production simulation (step 4). During these runs, the first 2.5 ns was discarded as additional equilibration, and the average and uncertainty (or error) was computed using the GROMACS tool g_energy71 by breaking the remaining 12.5 ns into five blocks of 2.5 ns each. For each species, simulations were performed at 1 bar and a range of temperatures. Temperatures were chosen that were expected to be below the normal boiling point and above the normal melting point. Estimates of these properties may be taken from experiment or prior molecular simulation results when available. Although not necessary in the present study, when unavailable, these properties may be estimated using the group contribution method of Hukkerikar et al.7 3.2.1. Free Energy Calculations. The residual chemical potential (μres) of the systems was computed using a multistage free energy perturbation approach19−21,32,82 with the multistate Bennett’s acceptance ratio method (MBAR).83−86 The main idea is to connect two systems of interest via a series of intermediate stages. The two systems of interest in the present study are a pure liquid in which all N molecules are fully interacting and the same pure liquid with N − 1 fully interacting molecules and a single molecule decoupled from (or noninteracting with) the system at the same temperature and pressure; to facilitate discussion, we will refer to the molecule that is being coupled/decoupled in our pure liquid as the solute. The change in Gibbs free energy between these two systems gives the desired residual chemical potential. These systems are connected through a series of stages that begin with a noninteracting solute molecule in the pure liquid (i.e., an ideal gas reference state at T and p′) and end with a fully interacting solute in the pure liquid. The intermediate stages serve to scale the intermolecular interaction potential of the solute. A specific stage is designated by index m. Intermolecular LJ and electrostatic interactions are regulated by the stage dependent elec coupling parameter λLJ m and λm , respectively, which vary from 0 LJ elec elec ≤ λm ≤ 1 and 0 ≤ λm ≤ 1. When λLJ m = λm = 1, the solute is LJ elec fully coupled to the system. When λm = λm = 0, the solute is decoupled from the system. To prevent instabilities in the trajectory when the solute molecule is nearly decoupled from D


Article



insertions/deletions.40 The simulations employed the coupleddecoupled CBMC formulation40,42,95,96 with parameters (or settings) from Martin and Frischknecht.97 Throughout every CBMC move (i.e., for all trial positions), the full nonbonded cutoff (rcut) was used. For each T, four independent simulations were performed. Each run was initiated with an empty box and with a different seed for the random number generator. The simulations were then analyzed independently. The reported results are the average value for the four runs, and the reported uncertainty corresponds to the 90% confidence interval.98 All properties were computed directly from the molecule number probability distribution. The simulations were run for between 25 and 325 × 106 MC steps as necessary to converge the molecule number probability distribution, with the TMMC weighting function updated every 1 × 106 MC steps.

4. RESULTS AND DISCUSSION 4.1. Alkanes and Alcohols. The goal of this work is to demonstrate the use of established molecular simulation free energy methods to compute the pure liquid fugacity of low volatile liquids. We therefore begin by examining n-octane, 3,4dimethylhexane, and cyclohexane all modeled using the TraPPE-UA force field. These systems were chosen because they are united atom models containing no partial charges allowing reference GC-TMMC results to readily be generated, and they have all been studied previously. While the intermolecular interaction potential may be deemed relatively simple, branching in 3,4-dimethylhexane and cyclic cyclohexane present challenges for CBMC methodologies.39,40,42,44,99 The liquid molar volume and fugacity are examined in Figures 2 and 3, respectively. For n-octane and cyclohexane, the MD and GC-TMMC results of the present study are in excellent agreement. The MD simulations were all performed at temperatures below the normal boiling point where we expect that f L ≈ fsat ≈ psat. In all cases, the GC-TMMC results show that as we approach the normal boiling point from above, fsat and psat become indistinguishable as expected, which in turn agrees with f L computed in our MD simulations. GC-TMMC calculations were performed down to 240 K (vapor pressure of 72 ± 3 Pa) for n-octane and down to 323 K (vapor pressure of 77 ± 2 kPa) for cyclohexane. Below these temperatures it was difficult to converge the molecule number probability distribution in a reasonable amount of simulation time because of the large free energy barriers separating the liquid and vapor phases. On the other hand, MD simulations were carried out at lower temperatures without a challenge. Note that recent MC methods have been developed to overcome this challenge (see for example refs 100 and 101) but are beyond the scope of the present study. While it is sufficient to show that the MD and GC-TMMC results of the present study are consistent for n-octane and cyclohexane, comparison is additionally made to reference simulation results.35,41,100 For all cases, the reference simulation results are in excellent agreement with the present study. Additionally, the low temperature MD results for n-octane and cyclohexane were found to be in remarkably good agreement with the NPT-TEE results of Rane et al.100 over the entire temperature range studied. The agreement between the MD and GC-TMMC results for 3,4-dimethylhexane was also good. The agreement between the MD and GC-TMMC results for the liquid fugacity is very good, which likewise shows agreement with the vapor pressure data of

Figure 2. Molar volume of n-octane, 3,4-dimethylhexane, and cyclohexane computed using molecular simulation with the TraPPEUA force field. Circles (○) and solid line with crosses (×) are MD and GC-TMMC results, respectively, from the present study. Triangles up (△) are Gibbs ensemble results from the original TraPPE-UA publications.41,42 Triangles down (▽) are GC-TMMC results from ref 35, and the broken lines are NpT-TEE results from ref 100. Note that the MD results are for the liquid phase at 1 bar while the others correspond to the liquid phase at saturation. The error in the data from the present study is smaller than the symbols and may be found in the Supporting Information of this manuscript.

Martin and Siepmann41 computed using MC calculations in the Gibbs ensemble. The liquid molar volumes, however, show slight deviations. The MD results are in reasonable agreement with the GC-TMMC results, although the MD results are slightly smaller. Likewise, at higher temperatures the GCTMMC results appear to be slightly smaller than Martin and Siepmann.41 As a result of the excellent agreement shown for noctane and cyclohexane, we suspect this slight discrepancy of the liquid molar volume is a consequence of the challenge associated with sampling molecules with branch points via Monte Carlo methods. Having demonstrating agreement between the free energy based approach (using MD in the present study) and standard MC methods for the studied alkanes, we next considered methanol, 1-propanol, and 2-propanol again modeled with the TraPPE-UA force field. As compared to the alkanes, the alcohols all contained partial charges (electrostatic interactions) and are known to have strong directional intermolecular interactions and structure in the liquid phase. Additionally, 2propanol contains a branch point. As shown in Figures 4 and 5, for all cases, the MD results of the present study and reference MC results are in excellent agreement.35,43,100 For the case of 1-propanol, the reference NPT-TEE results of Rane et al.100 span the same low temperature region considered in the present study. Once again, we find that the results are in remarkable agreement, down to a predicted vapor pressure of 3.7 ± 0.3 Pa at 225 K. This low vapor pressure would correspond to a vapor box with a density of approximately 93.40 nm3 per molecule. 4.2. Acetanilide, Acetaminophen, and Phenacetin. The present study is motivated by the desire to compute the activity coefficient of solid solutes in solution at ambient E


Article



Figure 3. Liquid fugacity and vapor pressure of n-octane, 3,4dimethylhexane, and cyclohexane computed using molecular simulation with the TraPPE-UA force field. Circles (○) and solid lines are MD and GC-TMMC results, respectively, from the present study. The MD results are for the pure liquid fugacity at 1 bar. The solid red line is GC-TMMC results for the pure fugacity at saturation, and the solid black line is GC-TMMC results for the vapor pressure at saturation. The reference data are all for the vapor pressure at saturation: triangles up (△) are Gibbs ensemble results from the original TraPPE-UA publications;41,42 triangles down (▽) are GC-TMMC results from ref 35; the broken lines are NpT-TEE results from ref 100. The dotted horizontal line corresponds to 1 atm and is drawn as a reference. The error in the data from the present study is smaller than the symbols and may be found in the Supporting Information of this manuscript.

Figure 5. Liquid fugacity and vapor pressure of methanol, 1-propanol, and 2-propanol computed using molecular simulation with the TraPPE-UA force field. Circles (○) correspond to the pure liquid fugacity at 1 bar computed using MD in this study. The reference data are all for the vapor pressure at saturation: triangles up (△) are Gibbs ensemble results from the original TraPPE-UA publications;43 triangles down (▽) are GC-TMMC results from ref 35; the broken line is NpT-TEE results from ref 100. The dotted horizontal line corresponds to 1 atm and is drawn as a reference. The error in the data from the present study is smaller than the symbols and may be found in the Supporting Information of this manuscript.

conditions. This therefore demands the calculation of the pure fugacity of a subcooled liquid. It is advisible not to model such systems directly, as they are quasi-nonergodic, requiring long simulation times and/or sophisticated enhanced sampling techniques to obtain accurate thermodynamic properties.29−32 Results for acetanilide, acetaminophen, and phenacetin are shown in Figures 6 and 7. All three compounds are solid at ambient conditions. Our initial simulations were therefore restricted to temperatures just below the experimental normal boiling point14,15,102,103 and near the experimental normal melting point.104 Experimental liquid molar volumes over this temperature range were difficult to find. A limited amount of data were found for acetanilide at temperatures just above the melting point.105,106 These experimental data are in excellent agreement with the MD results. This is very promising given the predictive nature of the MD results with the force field assembled in the present study. We next turn our attention to the computed pure liquid fugacities. Vapor pressure data for these compounds are more abundant than liquid molar volumes. For both acetanilide and phenacetin, the MD results are in relatively close agreement with the experimental results at the higher temperatures. The deviation between the two increases as the temperature decreases. This trend suggests that the present study overpredicts the enthalpy of vaporization relative to experiment. Nonetheless, given the purely predictive nature of the MD simulations, the agreement is quite good. For the case of acetanilide at 400 K, the MD results predict a pure liquid fugacity of approximately 9.4 ± 0.8 Pa while the experimental Antoine equation evaluated at this temperature gives a vapor pressure of 247 Pa.14,103 Recalling 1 atm = 101. 325 × 105 Pa and given the challenge of the experiments at these conditions,

Figure 4. Molar volume of methanol, 1-propanol, and 2-propanol computed using molecular simulation with the TraPPE-UA force field. Circles (○) are MD results from the present study. Triangles up (△) are Gibbs ensemble results from the original TraPPE-UA publications.43 Triangles down (▽) are GC-TMMC results from ref 35, and the broken line is NpT-TEE results from ref 100. Note that the MD results are for the liquid phase at 1 bar while the others correspond to the liquid phase at saturation. The error in the data from the present study is smaller than the symbols and may be found in the Supporting Information of this manuscript.

F


Article



Figure 7. Liquid fugacity and vapor pressure of acetanilide, phenacetin, and acetaminophen. Circles (○) correspond to the pure liquid fugacity at 1 bar computed using MD in this study. The solid black line is a Clausius−Clapeyron equation (ln f L = aT−1 + b) fit to the MD results. For acetaminophen, the fit is performed over the temperature range of 450−700 K. All other data are experimental data for the vapor pressure at saturation: triangles up (△) are from ref 102; solid red lines are the Antoine equations from ref 103; the solid green line is the Antoine equation for acetanilide from ref 14; the broken blue and green lines correspond to the Clausius−Clapeyron equation fit to liquid thermogravimetry data under nonisothermal and isothermal conditions from ref 15, respectively. The dotted horizontal line corresponds to 1 atm and the vertical dotted line corresponds to the normal melting point104 and are drawn as a reference. The error in the data from the present study is smaller than the symbols and may be found in the Supporting Information of this manuscript.

Figure 6. Molar volume of acetanilide, phenacetin, and acetaminophen at 1 bar. Circles (○) are MD results from the present study. Triangles up (△) and triangles down (▽) are experimental data from refs 105 and 106, respectively. The solid black line is a linear regression of the form ln v = aT + b fit to the MD results, where a and b are constants. For acetaminophen, the fit is performed over the temperature range of 450−700 K. The error in the data from the present study is smaller than the symbols and may be found in the Supporting Information of this manuscript.

the agreement is quite satisfying. As shown in the bottom pane of Figure 7, the computed liquid fugacity for acetaminophen agrees well with available experimental data,15 with the MD results slightly underpredicting the experimental data. To estimate the pure fugacity of subcooled acetanilide, phenacetin, and acetaminophen at ambient conditions, the MD results from the temperature range just below the experimental normal boiling point to near the experimental normal melting point were fit to a Clausius−Clapeyron equation of the form ln f L = aT−1 + b, where a and b are constants and the inverse of the uncertainty of f L was used to weight each datum during the regression. Below the normal boiling point, the Clausius− Clapeyron equation is expected to fit the data well.1 In all cases we obtain a coefficient of determination (R2 value) close to unity. The resulting expression could then be used to extrapolate the MD results to ambient conditions. Before changing topics, while it is satisfying that the MD results are in relatively good agreement with available experimental data, this is not the objective of the present study. More important is that we can compute the pure liquid fugacity above the normal melting point and extrapolate to ambient conditions. When coupled with free energy calculations for the solute of interest in solution, the activity coefficient may be computed in a self-consistent fashion. 4.2.1. Subcooled Liquids. If one desires the pure fugacity of subcooled liquid acetanilide, phenacetin, or acetaminophen at ambient conditions in order to compute an activity coefficient, it is tempting to directly perform the necessary free energy calculations at the temperature of interest. After all, this would result in a significant reduction of CPU and user time requiring just a single simulation. Acetaminophen is a current, popular pharmaceutical compound for which the present study is motivated. We therefore decided to perform MD simulations at

temperatures below the experimental normal melting point, with the lowest temperature being 300 K. The simulations were all performed following the protocol used for the higher temperatures consider. The results of these calculation are shown in Figure 7. We find that the direct simulation results deviate significantly from the Clausius−Clapeyron fit, with the deviation increasing as temperature decreases. In fact, we see that if these points were likewise fit to a Clausius−Clapeyron equation, we would find that the slope would differ from that above the experimental normal melting point. This would suggest a step change in the enthalpy of vaporization and the presence of a first order phase transition.107 Notice in Figure 6, however, that the molar volume remains a continuous function. Nucleation is a rare event; observing a liquid−solid phase transition during the finite simulation time is unrealistic. Nonetheless, the results of the present study are interesting and suggest the use of caution when performing similar calculations. We emphasize that the direct simulation results below the experimental melting point differ significantly from the Clausius−Clapeyron equation fit to data between the normal boiling point and experimental melting point extrapolated to ambient conditions. The Clausius−Clapeyron equation is expected to perform well for these conditions,1 suggesting the inappropriateness of the direct simulations. At 300 K, the direct simulation results and the extrapolated Clausius−Clapeyron equation differ by over 10 log units (or 10 kBT), which would cause significant differences in computed activity coefficients. G



■

The dynamics of a MD simulation are known to be dependent on the friction constant of the stochastic thermostat108 which was used in the present study to regulate the temperature. To ensure that this was not causing an issue here, the three simulations below the normal melting point were repeated with the friction constant for the stochastic thermostat decreased from 1.0 to 0.125 ps−1. This, however, had only a minute effect on the results. All of the molecular simulation results are tabulated in the Supporting Information of this manuscript.

ASSOCIATED CONTENT

* Supporting Information S

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.iecr.5b01827. Alternative derivation of the expression for the vapor pressure, the relationship between the Hildebrand solubility parameter and the computed fugacity, additional details about the GC-TMMC calculations, and the tabulated simulation results from the present study (PDF) All of the GROMACS force field files for acetanilide, acetaminophen, and phenacetin, catenated in ie5b01827_si_002.txt (TXT) All of the Gaussian 09 input files for acetanilide, acetaminophen, and phenacetin, catenated in ie5b01827_si_003.txt (TXT)

5. SUMMARY AND CONCLUSIONS


Article

In the current study, conventional molecular simulation free energy calculations and standard thermodynamic relations were applied to compute the pure liquid fugacity of low volatile liquids. The pure liquid fugacity is the normalization constant necessary to compute activity coefficients with respect to a Lewis−Randall (or Raoult’s Law) standard state. The ability to compute activity coefficients via molecular simulation may lead to the development of new simulation methodologies to predict solid−liquid equilibrium, in addition to allowing new solution theories and excess Gibbs free energy models to be developed. The method used in the present study is particularly attractive in that it uses the same molecular simulation free energy methods used to compute solvation free energies, which are available in many molecular dynamics and Monte Carlo software. This simplicity will allow the method to be used by novice practitioners, in addition to ensuring that activity coefficients are computed in a self-consistent manner. Calculations were performed for the reference compounds noctane, 3,4-dimethylhexane, cyclohexane, methanol, 1-propanol, and 2-propanol and excellent agreement was found between the presented method and reference Monte Carlo simulations. Calculations were additionally performed for acetanilide, phenacetin, and acetaminophen, all of which are solid at ambient conditions. These systems are of most interest and are also most challenging because computing the activity coefficient at ambient conditions requires knowledge of the fugacity of a subcooled liquid. Simulations of a subcooled liquid are not recommended as the system is quasi-nonergodic. We demonstrated how simulations may be performed at elevated temperatures and extrapolated to subcooled conditions. The results were encouraging, especially given the challenge of experimentally measuring the vapor pressure of compounds just above their normal melting point. While simulations involving a subcooled liquid are not advisable, it is tempting to do so. This is especially true for a novice practitioner for whom this work is targeted. It is the case that one may always perform a molecular simulation and one will always get some result. The question then is the truthfulness of the result. Here we performed calculations for acetaminophen below the experimental normal melting point, down to a temperature of 300 K. We found that the computed liquid fugacity deviated significantly from that extrapolated from higher temperatures. The value at 300 K differed by over 10 log units (or 10 kBT), which would cause significant differences in computed activity coefficients and ultimately computed compositions at equilibrium.

■

AUTHOR INFORMATION

Corresponding Author

*Phone: (513) 529-0784. Fax: (513) 529-0761. E-mail: [email protected]. Author Contributions †

G.B.F. and R.T.L. contributed equally to this work.

Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS A.S.P. gratefully acknowledges funding from the Miami University Committee on Faculty Research and start-up support from the Miami University College of Engineering and Computing. G.B.F. acknowledges summer support from the Miami University College of Engineering and Computing. Computing support was provided by the Ohio Supercomputer Center and Miami University’s Research Computing Support group. The authors also thank David L. Mobley, Pavel V. Klimovich, and Shuai Liu (University of California, Irvine) for assistance with PyMBAR.

■

REFERENCES

(1) Poling, B. E.; Prausnitz, J. M.; O’Connell, J. P. The Properties of Gases and Liquids, 5th ed.; McGraw-Hill: New York, NY, 2001. (2) Gani, R.; Muro-Suñ e,́ N.; Sales-Cruz, M.; Leibovici, C.; O’Connell, J. P. Mathematical and numerical analysis of classes of property models. Fluid Phase Equilib. 2006, 250, 1. (3) Gmehling, J. Present status and potential of group contribution methods for process development. J. Chem. Thermodyn. 2009, 41, 731. (4) O’Connell, J. P.; Gani, R.; Mathias, P. M.; Maurer, G.; Olson, J. D.; Crafts, P. A. Thermodynamic Property Modeling for Chemical Process and Product Engineering: Some Perspectives. Ind. Eng. Chem. Res. 2009, 48, 4619. (5) Harini, M.; Adhikari, J. A Review on Property Estimation Methods and Computational Schemes for Rational Solvent Design: A Focus on Pharmaceuticals. Ind. Eng. Chem. Res. 2013, 52, 6869. (6) Xue, Z.; Mu, T.; Gmehling, J. Comparison of the a Priori COSMO-RS Models and Group Contribution Methods: Original UNIFAC, Modified UNIFAC(Do), and Modified UNIFAC(Do) Consortium. Ind. Eng. Chem. Res. 2012, 51, 11809. (7) Hukkerikar, A. S.; Sarup, B.; Kate, A. T.; Abildskov, J.; Sin, G.; Gani, R. Group-contribution+ (GC+) based estimation of properties of pure components: Improved property estimation and uncertainty analysis. Fluid Phase Equilib. 2012, 321, 25. (8) Prausnitz, J. M.; Lichtenthaler, R. N.; de Azevedo, E. G. Molecular Thermodynamics of Fluid-Phase Equilibria, 2nd ed.; Prentice-Hall, Inc.: Englewood Cliffs, NJ, 1986. H


Article


Industrial & Engineering Chemistry Research (9) Hildebrand, J. H.; Prausnitz, J. M.; Scott, R. L. Regular and Related Solutions; Van Nostrand Reinhold Company: New York, NY, 1970. (10) Goodman, M. A. Vapor Pressure of Agrochemicals by the Knudsen Effusion Method Using a Quartz Crystal Microbalance. J. Chem. Eng. Data 1997, 42, 1227. (11) Monte, M. J. S.; Santos, L. M. N. B. F.; Fulem, M.; Fonseca, J. M. S.; Sousa, C. A. D. New Statis Apparatus and Vapor Pressure of Reference Materials: Naphthalene, Benzoic Acid, Benzophenone, and Ferrocene. J. Chem. Eng. Data 2006, 51, 757. (12) Goldfarb, J. L.; Suuberg, E. M. Vapor Pressures and Enthalpies of Sublimation of Ten Polycyclic Aromatic Hydrocarbonds Determined via Knudsen Effusion Method. J. Chem. Eng. Data 2008, 53, 670. (13) Booth, A. M.; Markus, T.; McFiggans, G.; Percival, C. J.; Mcgillen, M.; Topping, D. O. Design and construction of a simple Knudsen Effusion Mass Spectrometer (KEMS) system for vapor pressure measurements of low volatility organics. Atmos. Meas. Tech. 2009, 2, 355. (14) Vecchio, S.; Catalani, A.; Rossi, V.; Tomassetti, M. Thermal analysis study on vaporization of some analgesics. Acetanilide and derivatives. Thermochim. Acta 2004, 420, 99. (15) Vecchio, S.; Tomassetti, M. Vapor pressures and standard molar enthalpies, entropies and Gibbs energies of sublimation of three 4substituted acetanilide derivatives. Fluid Phase Equilib. 2009, 279, 64. (16) Farhod Chasib, K. Extraction of Phenolic Pollutants (Phenol and p-Chlorophenol) from Industrial Wastewater. J. Chem. Eng. Data 2013, 58, 1549. (17) Lapworth, D. J.; Baran, N.; Stuart, M. E.; Ward, R. S. Emerging organic contaminents in groundwater: A review of sources, fate and occurance. Environ. Pollut. 2012, 163, 287. (18) Gao, H.; Li, Y.; Wu, Y.; Luo, M.; Li, Q.; Xing, J.; Liu, H. Extractive Desulfurization of Fuel Using 3-Methylpyridinium-Based Ionic Liquids. Energy Fuels 2009, 23, 2690. (19) Shing, K. S.; Chung, S. T. Computer Simulation Methods for the Calculation of Solubility in Supercritical Extraction Systems. J. Phys. Chem. 1987, 91, 1674. (20) Kofke, D. A.; Cummings, P. T. Quantitative comparison and optimization of methods for evaluating the chemical potential by molecular simulation. Mol. Phys. 1997, 92, 973. (21) Shirts, M. R.; Pitera, J. W.; Swope, W. C.; Pande, V. S. Extremely precise free energy calculations of amino acid side chain analogs: Comparison of common molecular mechanics force fields for proteins. J. Chem. Phys. 2003, 119, 5740. (22) Lin, S. T.; Chang, J.; Wang, S.; Goddard, W. A., III; Sandler, S. I. Prediction of Vapor Pressures and Enthalpies of Vaporization Using a COSMO Solvation Model. J. Phys. Chem. A 2004, 108, 7429. (23) Winget, P.; Hawkins, G. D.; Cramer, C. J.; Truhlar, D. G. Predicting the Vapor Pressures from Self-Solvation Free Energies Calculated by the SM5 Series of Universal Solvation Models. J. Phys. Chem. B 2000, 104, 4726. (24) Möller, D.; Fischer, J. Vapour liquid equilibrium of a pure fluid from test particle method in combination with NpT molecular dynamics simulations. Mol. Phys. 1990, 69, 463. (25) Lotfi, A.; Vrabec, J.; Fischer, J. Vapour liquid equilibria of the Lennard-Jones fluid from the NpT plus test particle method. Mol. Phys. 1992, 76, 1319. (26) Panagiotopoulos, A. Z. Monte Carlo methods for phase equilibria of fluids. J. Phys.: Condens. Matter 2000, 12, R25. (27) Horn, H. W.; Swope, W. C.; Pitera, J. W. Characterization of the TIP4P-Ew water model: Vapor pressure and boiling point. J. Chem. Phys. 2005, 123, 194504. (28) Ahmed, A.; Sandler, S. I. Physicochemical Properties of Hazardous Energetic Compounds from Molecular Simulation. J. Chem. Theory Comput. 2013, 9, 2389. (29) Barrat, J. L.; Klein, M. L. Molecular Dynamics Simulations of Supercooled Liquids Near the Glass Transition. Annu. Rev. Phys. Chem. 1991, 42, 23.

(30) Mountain, R. D.; Thirumalai, D. Loss of ergodicity in glassy systems. AIP Conf. Proc. 1991, 256, 165. (31) Thirumalai, D.; Mountain, R. D.; Kirkpatrick, T. R. Ergodic behavior in supercooled liquids and in glasses. Phys. Rev. A: At., Mol., Opt. Phys. 1989, 39, 3563. (32) Chipot, C., Pohorille, A., Eds. Free Energy Calculations: Theory and Applications in Chemistry and Biology; Springer Series in Chemical Physics, Vol. 86; Springer: New York, NY, 2007. (33) Errington, J. R. Direct calculation of liquid-vapor phase equilibria from transition matrix Monte Carlo. J. Chem. Phys. 2003, 118, 9915. (34) Singh, J. K.; Errington, J. R. Calculation of Phase Coexistence Properties and Surface Tensions of n-Alkanes with Grand-Canonical Transition-Matrix Monte Carlo Simulation and Finite-Size Scaling. J. Phys. Chem. B 2006, 110, 1369. (35) Paluch, A. S.; Shen, V. K.; Errington, J. R. Comparing the Use of Gibbs Ensemble and Grand-Canonical Transition-Matrix Monte Carlo Methods to Determin Phase Equilibria. Ind. Eng. Chem. Res. 2008, 47, 4533. (36) Ferrenberg, A. M.; Swendsen, R. H. New Monte Carlo Technique for Studying Phase Transitions. Phys. Rev. Lett. 1988, 61, 2635. (37) Errington, J. R.; Panagiotopoulos, A. Z. Phase equilibria of the modified Buckingham exponential-6 potential from Hamiltonian scaling grand canonical Monte Carlo. J. Chem. Phys. 1998, 109, 1093. (38) Hill, T. L. Statistical Mechanics: Principles and Selected Applications; Dover Publications, Inc.: Mineola, NY, 1987. (39) Macedonia, M. D.; Maginn, E. J. A biased grand canonical Monte Carlo method for simulating adsorption using all-atom and branched united atom models. Mol. Phys. 1999, 96, 1375. (40) Frenkel, D.; Smit, B. Understanding Molecular Simulation: From Algorithms to Applications, 2nd ed.; Academic Press: San Diego, CA, 2002. (41) Martin, M. G.; Siepmann, J. I. Transferable Potentials for Phase Equilibria. 1. United-Atom Description of n-Alkanes. J. Phys. Chem. B 1998, 102, 2569. (42) Martin, M. G.; Siepmann, J. I. Novel Configurational-Bias Monte Carlo Method for Branced Molecules. Transferable Potentials for Phase Equilibria. 2. United-Atom Description of Branched Alkanes. J. Phys. Chem. B 1999, 103, 4508. (43) Chen, B.; Potoff, J. J.; Siepmann, J. I. Monte Carlo Calculations for Alcohols and Their Mixtures with Alkanes. Transferable Potentials for Phase Equilibria. 5. United-Atom Description of Primary, Secondary, and Tertiary Alcohols. J. Phys. Chem. B 2001, 105, 3093. (44) Lee, J. S.; Wick, C. D.; Stubbs, J. M.; Siepmann, J. I. Simulating the vapour-liquid equilibria of large cyclic alkanes. Mol. Phys. 2005, 103, 99. (45) Rai, N.; Siepmann, J. I. Transferable Potentials for Phase Equilibria. 9. Explicit Hydrogen Description of Benzene and FiveMembered and Six-Membered Heterocyclic Aromatic Compounds. J. Phys. Chem. B 2007, 111, 10790. (46) Rai, N.; Bhatt, D.; Siepmann, J. I.; Fried, L. E. Monte Carlo simulations of 1,3,5-triamino-2,4,6-trinitrobenzene (TATB): Pressure and temperature effects for the solid phase and vapor-liquid phase equilibria. J. Chem. Phys. 2008, 129, 194510. (47) Rai, N.; Siepmann, J. I. Transferable Potentials for Phase Equilibria. 10. Explicit-Hydrogen Description of Substituted Benzenes and Polycyclic Aromatic Compounds. J. Phys. Chem. B 2013, 117, 273. (48) Wick, C. D.; Stubbs, J. M.; Rai, N.; Siepmann, J. I. Transferable Potentials for Phase Equilibria. 7. Primary, Secondary, and Tertiary Amines, Nitroalkanes and Nitrobenzene, Nitriles, Amides, Pyridine, and Pyrimidine. J. Phys. Chem. B 2005, 109, 18974. (49) Note that TraPPE-EH does parametrize aniline in ref 46. We chose not to use the N LJ parameters from aniline because it is a primary amine, and it was previously shown in ref 48 that there is an appreciable change in LJ parameters in going from a primary to secondary amine. Also, in ref 48, the primary amide LJ parameters are the same as the primary amine (i.e., only the charges change). Therefore, we took the secondary amide LJ parameters for N to be the I


Article


Industrial & Engineering Chemistry Research same as for a secondary amine. The H has LJ parameters of 0 in primary and secondary amines and primary amides. The LJ parameters for H in a secondary amide were therefore taken to also be 0 in the present study. (50) Stubbs, J. M.; Potoff, J. J.; Siepmann, J. I. Transferable Potentials for Phase Equilibria. 6. United-Atom Description for Ethers, Glycols, Ketones, and Aldehydes. J. Phys. Chem. B 2004, 108, 17596. (51) The alkyl groups (CH3 and CH2) were modeled as a single united-atom pseudoatom as a result of the parametrization of the TraPPE-EH force field for n-alkanes, which places the LJ site for a hydrogen atom at the center of the corresponding bond,109 and the complication of implementing such a model in a molecular dynamics framework. (52) Zhao, Y.; Truhlar, D. G. The M06 theory of density functionals for main group thermochemistry, thermochemical kinetics, noncovalent interactions, excited states, and transition elements: two new functionals and systematic testing of four M06-class functionals and 12 other functionals. Theor. Chem. Acc. 2008, 120, 215. (53) Cramer, C. J. Essentials of Computational Chemistry; John Wiley & Sons Ltd: Chichester, West Sussex, England, 2002. (54) Marenich, A. V.; Cramer, C. J.; Truhlar, D. G. Universal solvation model based on solute electron density and on a continuum model of the solvent defined by the bulk dielectric constant and atomic surface tensions. J. Phys. Chem. B 2009, 113, 6378. (55) Marenich, A. V.; Cramer, C. J.; Truhlar, D. G. Performance of SM6, SM8, and SMD on the SAMPL1 test set for the prediction of small-molecule solvation free energies. J. Phys. Chem. B 2009, 113, 4538. (56) Ribeiro, R. F.; Marenich, A. V.; Cramer, C. J.; Truhlar, D. G. Prediction of SAMPL2 aqueous solvation free energies and tautomeric ratios using the SM8, SM8AD, and SMD solvation models. J. Comput.Aided Mol. Des. 2010, 24, 317. (57) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, J. A.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, J. M.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, Ö ; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Gaussian 09, revision A.02; Gaussian, Inc.: Wallingford, CT, 2009. (58) Bayly, C. I.; Cieplak, P.; Cornell, W. D.; Kollman, P. A. A WellBehaved Electrostatic Potential Based Method Using Charge Restraints for Deriving Atomic Charges: The RESP Model. J. Phys. Chem. 1993, 97, 10269. (59) Cieplak, P.; Cornell, W. D.; Bayly, C.; Kollman, P. A. Application of the multimolecule and multiconformational RESP methodology to biopolymers: Charge derivation for DNA, RNA, and proteins. J. Comput. Chem. 1995, 16, 1357. (60) Case, D. A.; Cheatham, T.; Darden, T.; Gohlke, H.; Luo, R.; Merz, K. M.; Onufriev, A.; Simmerling, C.; Wang, B.; Woods, R. The Amber biomolecular simulation programs. J. Comput. Chem. 2005, 26, 1668. (61) Case, D. A.; Darden, T. A.; Cheatham, T. E., III; Simmerling, C. L.; Wang, J.; Duke, R. E.; Luo, R.; Walker, R. C.; Zhang, W.; Merz, K. M.; Roberts, B.; Hayik, S.; Roitberg, A.; Seabra, G.; Swails, J.; Götz, A. W.; Kolossváry, I.; Wong, K. F.; Paesani, F.; Vanicek, J.; Wolf, R. M.; Liu, J.; Wu, X.; Brozell, S. R.; Steinbrecher, T.; Gohlke, H.; Cai, Q.; Ye, X.; Wang, J.; Hsieh, M. J.; Cui, G.; Roe, D. R.; Mathews, D. H.; Seetin, M. G.; Salomon-Ferrer, R.; Sagui, C.; Babin, V.; Luchko, T.; Gusarov, S.; Kovalenko, A.; Kollman, P. AMBER 12; University of California San Francisco: San Francisco, CA, 2012.

(62) Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. Development and Testing of a General Amber Force Field. J. Comput. Chem. 2004, 25, 1157. (63) Wang, J.; Wang, W.; Kollman, P. A.; Case, D. A. Automatic Atom Type and Bond Type Perception in Molecular Mechanical Calculations. J. Mol. Graphics Modell. 2006, 25, 247. (64) Sousa da Silva, A. W.; Vranken, W. F. ACPYPE - AnteChamber PYthon Parser interfacE. BMC Res. Notes 2012, 5, 367. (65) Sousa da Silva, A. W.; Vranken, W. F. ACPYPE: AnteChamber PYthon Parser interfacE. https://code.google.com/p/acpype/ (accessed May 1, 2014). (66) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 435. (67) Pronk, S.; Páll, S.; Schulz, R.; Larsson, P.; Bjelkmar, P.; Apostolov, R.; Shirts, M. R.; Smith, J. C.; Kasson, P. M.; van der Spoel, D.; Hess, B.; Lindahl, E. GROMACS 4.5: a high-throughput and highly parallel open source molecular simulation toolkit. Bioinformatics 2013, 29, 845. (68) GROMACS: Fast, Flexible, Free. http://www.gromacs.org/ (accessed Aug 1, 2013). (69) Martínez, L.; Andrade, R.; Birgin, E. G.; Martínez, J. M. Packmol: A package for building initial configurations for molecular dynamics simulations. J. Comput. Chem. 2009, 30, 2157. (70) Packmol: Packing Optimization for Molecular Dynamics Simulations. http://www.ime.unicamp.br/martinez/packmol/ (accessed Aug 1, 2013). (71) van der Spoel, D.; Lindahl, E.; Hess, B.; The GROMACS Development Team. GROMACS User Manual, version 4.6.3; University of Groningen: Groningen, The Netherlands, 2013; ftp:// ftp.gromacs.org/pub/manual/manual-4.6.3.pdf. (72) van Gunsteren, W. F.; Berendsen, H. J. C. A leap-frog algorithm for stochastic dynamics. Mol. Simul. 1988, 1, 173. (73) Berendsen, H. J. C. Simulating the Physical World: Hierarchial Modeling from Quantum Mechanics to Fluid Dynamics; Cambridge University Press: New York, NY, 2007. (74) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. Molecular dynamics with coupling to an external bath. J. Chem. Phys. 1984, 81, 3684. (75) Parrinello, M.; Rahman, A. Polymorphic transitions in single crystals: A new molecular dynamics method. J. Appl. Phys. 1981, 52, 7182. (76) Nosé, S.; Klein, M. L. Constant pressure molecular dynamics for molecular systems. Mol. Phys. 1983, 50, 1055. (77) Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. G. E. M. LINCS: A linear constraint solver for molecular simulations. J. Comput. Chem. 1997, 18, 1463. (78) Hess, B. P-LINCS: A Parallel Linear Constraint Solver for molecular simulation. J. Chem. Theory Comput. 2008, 4, 116. (79) Leach, A. R. Molecular Modelling: Principles and Applications, 2nd ed.; Pearson Education Limited: Harlow, England, 2001. (80) Deserno, M.; Holm, C. How to mesh up Ewald sums. I. A theoretical and numerical comparison of various particle mesh routines. J. Chem. Phys. 1998, 109, 7678. (81) The alkanes n-octane, 3,4-dimethylhexane, and cyclohexane were all modeled with the TraPPE-UA alkane force field and contain no partial charges. For these systems, use of SPME for long-range electrostatic interactions is not necessary. Additionally, since the system have no charges, stages to decouple the electrostatic interactions are not necessary. (82) Kofke, D. A.; Cummings, P. T. Precision and accuracy of staged free-energy perturbation methods for computing the chemical potential by molecular simulation. Fluid Phase Equilib. 1998, 150, 41. (83) Bennett, C. H. Efficient Estimation of Free Energy Differences from Monte Carlo Data. J. Comput. Phys. 1976, 22, 245. (84) Shirts, M. R.; Bair, E.; Hooker, G.; Pande, V. S. Equilibrium Free Energies from Nonequilibrium Measurements Using MaximumLikelihood Methods. Phys. Rev. Lett. 2003, 91, 140601. J


Article


Industrial & Engineering Chemistry Research (85) Lu, N.; Singh, J. K.; Kofke, D. A. Appropriate methods to combine forward and reverse free-energy perturbation averages. J. Chem. Phys. 2003, 118, 2977. (86) Shirts, M. R.; Chodera, J. D. Statistically optimal analysis of samples from multiple equilibrium states. J. Chem. Phys. 2008, 129, 124105. (87) Beutler, T. C.; Mark, A. E.; van Schaik, R. C.; Gerber, P. R.; van Gunsteren, W. F. Avoiding singularities and numerical instabilities in free energy calculations based on molecular simulations. Chem. Phys. Lett. 1994, 222, 529. (88) Shirts, M. R.; Pande, V. S. Solvation free energies of amino acid side chain analogs for common molecular mechanics water models. J. Chem. Phys. 2005, 122, 134508. (89) Steinbrecher, T.; Mobley, D. L.; Case, D. A. Nonlinear scaling schemes for Lennard-Jones interactions in free energy calculations. J. Chem. Phys. 2007, 127, 214108. (90) PyMBAR: Python implementation of the multistate Bennett acceptance ratio (MBAR). https://github.com/choderalab/pymbar (accessed May 1, 2014). (91) Chodera, J. D.; Swope, W. C.; Pitera, J. W.; Seok, C.; Dill, K. A. Use of the Weighted Histogram Analysis Method for the Analysis of Simulated and Parallel Tempering Simulations. J. Chem. Theory Comput. 2007, 3, 26. (92) Klimovich, P. V.; Shirts, M. R.; Mobley, D. L. Guidelines for analysis of free energy calculations. J. Comput.-Aided Mol. Des. 2015, 29, 397. (93) Martin, M. G. MCCCS Towhee: a tool for Monte Carlo molecular simulation. Mol. Simul. 2013, 39, 1212. (94) Martin, M. G. MCCCS Towhee. http://towhee.sourceforge.net (accessed May 1, 2014) . (95) Martin, M. G.; Thompson, A. P. Industrial property prediction using Towhee and LAMMPS. Fluid Phase Equilib. 2004, 217, 105. (96) Martin, M. G.; Biddy, M. J. Monte Carlo Molecular Simulation Predictions for the Heat of Vaporization of Acetone and Butyramide. Fluid Phase Equilib. 2005, 236, 53. (97) Martin, M. G.; Frischknecht, A. L. Using arbitrary trial distributions to improve intramolecular sampling in configurationalbias Monte Carlo. Mol. Phys. 2006, 104, 2439. (98) Harris, D. C. Quantitative Chemical Analysis, 6th ed.; W. H. Freeman and Company: New York, NY, 2002; pp 66−68. (99) Shah, J. K.; Maginn, E. J. A general and efficient Monte Carlo method for sampling intramolecular degrees of freedom of branched and cyclic molecules. J. Chem. Phys. 2011, 135, 134121. (100) Rane, K. S.; Murali, S.; Errington, J. R. Monte Carlo Simulation Methods for Computing Liquid-Vapor Saturation Properties of Model Systems. J. Chem. Theory Comput. 2013, 9, 2552. (101) Rane, K. S.; Errington, J. R. Using Monte Carlo Simulation to Compute Liquid-Vapor Saturation Properties of Ionic Liquids. J. Phys. Chem. B 2013, 26, 8018. (102) Stull, D. R. Vapor Pressure of Pure Substances- Organic Compounds. Ind. Eng. Chem. 1947, 39, 517. (103) Stephenson, R. M.; Malanowski, S. Handbook of the Thermodynamics of Organic Compounds; Elsevier Science Publishing: New York, NY, 1987. Acetanilide is on p 264, and phenacetin is on page 326. (104) Marrero, J., Abildskov, J., Eds. Solubility and Related Properties of Large Complex Chemicals Part 1: Organic Solutes Ranging from C4 to C40; DECHEMA: Frankfurt am Main, Germany, 2003. (105) Turner, W. E. S.; Merry, E. W. CCXVIII.- The molecular complexity, in the liquid state, of tervalent nitrogen compounds. J. Chem. Soc., Trans. 1910, 97, 2069. (106) Dunstan, A. E.; Mussell, A. G. CCV.- The viscosity of certain amides. J. Chem. Soc., Trans. 1910, 97, 1935. (107) Hansen, J. P.; McDonald, I. R. Theory of Simple Liquids; Academic Press, Inc.: New York, NY, 1976. (108) Basconi, J. E.; Shirts, M. R. Effects of Temperature Control Algorithms on Transport Properties and Kinetics in Molecular Dynamics Simulations. J. Chem. Theory Comput. 2013, 9, 2887.

(109) Chen, B.; Siepmann, J. I. Transferable Potentials for Phase Equilibria. 3. Explicit-Hydrogen Description of Normal Alkanes. J. Phys. Chem. B 1999, 103, 5370.

K


Calculating the Fugacity of Pure, Low Volatile Liquids via Molecular

Recommend Documents