Computational Study of Molecular Structure and Self-Association of Tri

Jan 7, 2015 - Junju Mu , Ryuhei Motokawa , Kazuhiro Akutsu , Shotaro Nishitsuji , and Andrew J. Masters. The Journal of Physical Chemistry B 2018 122 ...
0 downloads 0 Views 2MB Size
Download PDF
Article pubs.acs.org/JPCB

Computational Study of Molecular Structure and Self-Association of Tri‑n‑butyl Phosphates in n‑Dodecane Quynh N. Vo,† Cory A. Hawkins,† Liem X. Dang,‡ Mikael Nilsson,† and Hung D. Nguyen*,† †

Department of Chemical Engineering and Materials Science, University of California, Irvine, Irvine, California 92697, United States Physical Science Division, Pacific Northwest National Laboratory, Richland, Washington 93352, United States

‡

S Supporting Information *

ABSTRACT: Tri-n-butyl phosphate (TBP) is an important extractant used in the solvent extraction process for recovering uranium and plutonium from used nuclear fuel. An atomistic molecular dynamics study was used to understand the fundamental molecularlevel behavior of extracting agents in solution. Atomistic parametrization was carried out using the AMBER force field to model the TBP molecule and n-dodecane molecule, a commonly used organic solvent. Validation of the optimized force field was accomplished through various thermophysical properties of pure TBP and pure n-dodecane in the bulk liquid phase. The mass density, dipole moment, self-diffusion coefficient, and heat of vaporization were calculated from our simulations and compared favorably with experimental values. The molecular structure of TBPs in n-dodecane at a dilute TBP concentration was examined based on radial distribution functions. 1D and 2D potential mean force studies were carried out to establish the criteria for identifying TBP aggregates. The dimerization constant of TBP in the TBP/n-dodecane mixture was also obtained and matched the experimental value.

1. INTRODUCTION In the race for sustainable energy sources, nuclear power is a strong candidate because it operates on the naturally abundant uranium and has the capability to provide a round-the-clock base load of power with no CO2 emissions. Despite the high potential, management of the used nuclear fuel is crucial. An effective strategy in dealing with the used fuel is to recycle certain elements in the material using advanced chemical separation, which can significantly reduce the high-level waste volume.1,2 In a typical solvent extraction (SX) process for the recovery and purification of metal ions from used nuclear fuel, the aqueous phase, containing the metal ions, is vigorously mixed with extracting agents and hydrocarbon diluent of the organic phase. This is followed by a separation process, in which the organic phase is segregated away from the aqueous phase. The transfer of metal ions from the aqueous phase into the organic solvent by such means allows straightforward implementation for continuous processes such as the phase separation process by the centrifugal contactors. The most commonly used SX technique is called PUREX (Plutonium Uranium Reduction EXtraction), which employs tri-n-butyl phosphate (TBP) as an extracting agent to selectively complex uranium and plutonium from the constituents of the used fuel.3,4 However, the metal extraction mechanism and phase separation involve many molecular-level events that are not well understood. Specifically, TBP tends to self-assemble into reverse micelles, in which the polar hydrophilic head groups and solutes (e.g., waters, acids) are surrounded by a shell comprised of the nonpolar hydrophobic tail groups, possibly interfering with the phase behavior. To intensify the problem, © XXXX American Chemical Society

these reverse micelles can aggregate and form an undesired third phase, which can interrupt the counter current SX process operation and can be potentially hazardous from a criticality point of view as metal ions may be concentrated in a reduced volume.5−9 Molecular dynamics (MD) studies have been employed to gain insight and to visualize molecular-level behavior of extraction systems.10−17 However, to date, accurate force fields for modeling TBP and hydrocarbon molecules in the liquid phase are still lacking. Previously, TBP had been examined as a function of its environments, in the gas phase, in chloroform, and in pure water by Wipff and co-workers who carried out MD studies.10 They employed the MNDO semiempirical method to obtain the partial atomic charges and the AMBER91 force field for TBP molecule. This study did not focus on TBP bulk behavior in the liquid phase or in organic solvent, which is the condition of the PUREX process. A recent comparative study from Khomami and co-workers investigated four different atomistic force fields including the one from Wipff et al.10 to study TBP bulk liquid behavior.15 Their results showed deviations of the mass density and dipole moment of TBP from available experimental data unless charges were significantly reduced by up to 40% of the original values. The self-diffusion coefficient (SDC) of TBP was considerably underpredicted by all force fields even with scaled charges. Received: October 14, 2014 Revised: January 5, 2015

A

DOI: 10.1021/jp510365c J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B Hence, there has not been a definitive conclusion of an accurate force field to be employed for the TBP molecule. The PUREX process uses various kerosene-type solvents, which are most commonly emulated by n-dodecane, a molecule that has been parametrized for MD simulations using different force fields. As a pure diluent, n-dodecane is known as a liquid at room temperature; however, previous studies reported the liquid-to-gel-phase transition at room temperature for longchain hydrocarbons of more than eight carbons (including ndodecane) using the OPLS-AA force field.18,19 Khomami and co-workers attempted to correct this behavior by changing the scaling factor (SF) for the 1−4 intramolecular electrostatic and van der Waals interaction.19 The default SF for both interactions in the OPLS-AA force field is 0.5; therefore, they tested the SF values of 0.4, 0.3, 0.2, 0.1, and 0. The simulated melting point of n-dodecane from their study did not match the experimental value except for a SF of 0, where the 1−4 interaction was completely disregarded. The calculated heat of vaporization and density of n-dodecane still deviated slightly while SDC was remarkably smaller than experimental results for all values of SF. Nevertheless, they used a SF of 0.3 or 0.4 for their n-dodecane model in later papers.16,17 Moreover, Bockmann and co-workers extensively reparameterized the OPLS-AA force field for alkanes and alkenes by refining the depth of the Lennard-Jones potential, torsional parameters, and partial charges, which resulted in improved heats of vaporization values, densities, and phase transition temperatures.18 In this work, we carried out a rigorous, self-consistent, and systematic characterization study on the bulk liquid behavior of TBP and n-dodecane separately by refining the Lennard-Jones (LJ) potential parameters of the general AMBER force field (GAFF)20 in order to achieve more accurate models for both TBP and n-dodecane molecules. We optimized the LJ potential parameters based on the experimental mass density and heat of vaporization of each molecule. Other properties such as SDC of pure TBP and n-dodecane and the dipole moment of TBP molecules were calculated and compared with literature values as a means of validation for the optimized LJ potential parameters. Because there was only one value for the SDC of TBP molecules reported in the literature back in 1953,21,22 we verified it by carrying out our own experiment adopting an NMR technique. Computational prediction on the density of TBP as a function of temperature was also carried out for comparison where experimental data was available. The best models of TBP and n-dodecane were then employed for long MD simulations at a low concentration. Radial distribution functions (RDFs) and the Potential Mean Force (PMF) were constructed to gain insights into the conformations and association of the TBP molecule with itself and the n-dodecane diluent. This allowed us to calculate the self-association constant of TBP molecules in n-dodecane from MD simulations for comparison with experimental results.

Figure 1. Schematic representations of (A) TBP and (B) n-dodecane molecules.

All simulations were performed using the GPU version of PMEMD of Amber12.27,28 The initial configurations of the studied systems were generated using the Packmol software package,29 where all of the molecules were randomly distributed. Equations of motion were integrated with the velocity Verlet algorithm with a time step of 1 fs and a cutoff of 10 Å. Long-range electrostatics (>10 Å) was treated by the particle-mesh Ewald (PME) method. The nonbonded interactions calculated for pairs of atoms in TBP and n-dodecane separated by three bonds (1−4 interaction) were applied with Amber default SFs 1/2 and 1/1.2 for van der Waals and electrostatics, respectively. Bonds involving hydrogen atoms were constrained using the SHAKE algorithm with a tolerance of 10−7. Three-dimensional cubic periodic boundary conditions were applied. Translational center of mass motions of the molecules were removed every 1000 steps (1 ps). The initial systems were subjected to energy minimization using the steepest decent algorithm until the energies equilibrated. Each system then underwent a heating process in which its temperature was increased gradually from 0 to 400 K for 1 ns. For the study of pure TBP at different temperatures, those systems with temperatures higher than 298 K were heated up to 600 K over 2 ns. All systems were then cooled slowly from 400 K (or 600 K) to the desired temperatures in 0.5 ns (or 1 ns) and held there for another 5−15 ns depending on the studied physical properties.30 The isothermal−isobaric (NPT) ensemble coupled with the Langevin thermostat (ntt = 3) and Berendsen barostat was used to converge the systems’ densities, temperatures, and pressures. The canonical ensemble (NVT), where the temperature was controlled by the weak-coupling algorithm (ntt = 1), was used for production runs. In the case of SDC studies, the microcanonical ensemble (NVE) with Ewald coefficient of 0.35 for the PME method was applied. Gas-phase simulations were modeled by using a single molecule of interest with no box information. The same constraints and simulation method were applied as in the liquid phase. 2.1. Parameterization of TBP and n-Dodecane. The parametrizations of TBP and n-dodecane were initialized using the default GAFF set. The densities and enthalpies of vaporization were calculated for each molecule and compared with experimental data. The LJ potential parameters, namely, σ and ε, were then refined accordingly to reproduce these physical properties of both molecules using the same optimization technique as reported in a previous study31 but for a nonpolarizable force field. Other properties such as the electric dipole moment (EDM) and SDC were used as a means to validate the refined force fields for TBP and n-dodecane. Different studied LJ parameter sets of both molecules can be

2. COMPUTATIONAL DETAILS A schematic representation of the TBP and n-dodecane molecules is shown in Figure 1. The GAFF20 parameter set for bond length stretching, bond angle bending, dihedral angle torsion, and nonbonded interactions was used for both TBP and n-dodecane. The partial electric charges of atoms of TBP and n-dodecane molecules were obtained from the AM1-BCC model23,24 using the antechamber 20,25 package that is implemented in Amber12.26 B

DOI: 10.1021/jp510365c J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

where pi is the instantaneous EDM of the ith molecule, qa is the charge of atom a, rCOM is the center of mass of the ith molecule, i and N is the total number of molecules in the system. 2.4. SDC. The SDC was obtained from the slope of the mean-square displacement of the center of mass averaged over the trajectories of individual molecules using the Einstein relation

found in the Supporting Information (Tables S1 and S2 and Figures S1 and S2). 2.2. Enthalpy of Vaporization. The enthalpy of vaporization of TBP and n-dodecane was computed using the equation32 potential potential ΔH vap(T ) = Egas (T ) − E liquid (T ) + RT

N

where Epotential and Epotential are the potential energies in the gas gas liquid and liquid phases, respectively; R is the gas constant; and T is the simulation temperature. The experimental value of the enthalpy of vaporization of TBP is not directly available. Instead, the value was obtained by fitting the TBP vapor pressure over a range of temperatures using the Clausius−Clapeyron equation with the assumption that the enthalpy of vaporization remains constant over the investigated range of temperatures (15−30 °C). ln P = −

ΔH vap RT

Ds =

where r(t) is the position vector of atom i at time t, N is the total number of atoms, and ⟨⟩ denotes the ensemble average. A linear relation for the mean-squared displacement (MSD) profile can only be achieved after extended simulation durations. Provided that long enough trajectories were used, the MSD method had been studied to be more advantageous over the velocity autocorrelation function following the Green− Kubo formula.40 Yeh and Hummer had shown that Ds obtained from simulations under periodic boundary conditions suffered from system size effects due to the simplification of long-range interaction calculations.41 A corrected value of Dscorr could be obtained from extrapolating the slope fitted to the Ds as a function of inverse box lengths to the infinite system limit (1/L → 0). For these reasons, we studied three different system sizes (N = 200, 350, and 500 molecules) for pure n-dodecane and four different system sizes (N = 200, 350, 500, and 1000 molecules) for pure TBP. Each MD trajectory was 100 ns long for the TBP molecule and 55 ns long for the n-dodecane molecule with the first 15 ns using the NPT ensemble and followed by the NVE ensemble for the rest of the trajectory. Simulation Ds was compared to a previously reported value obtained by capillary experiments.21,22 Moreover, we conducted additional NMR experiments using 1H DOSY. 1H NMR measurements employed an ADVANCE600 Bruker NMR, equipped with a 5 mm triple-resonance broad-band probe with z-axis gradients and the temperature maintained at 298 K. 1D and pseudo-2D experiments were performed using a longitudinal encode−decode stimulated echo pulse sequence with bipolar gradient pulse pairs as developed by Wu et al.42 The linear gradient ramp consisted of 16 steps from 2 to 95% attenuation, and the gradient strength was externally calibrated using neat n-pentanol (≥99% Sigma-Aldrich, St. Louis, MO). Nondeuterated samples required manual shimming of the zcomponents without a solvent lock. To avoid resonance distortion due to signal saturation, neat samples of TBP (puriss, ≥99.0%, Fluka Analytical, Germany) and n-pentanol were prepared in 2 mm thin-wall tubes and housed in standard 5 mm thin-wall tubes using glass spacers (Wilmad, Vineland, NJ). Diffusion coefficients were calculated automatically by the TopSpin software (Bruker Corp., v. 3.2) by fitting the peak intensity to the exponential42

+C

We chose to compare our simulated heat of vaporization of TBP to the data set provided by Skene et al.33 They reported six different sets of vapor pressure for TBP over the temperature range of 0−40 °C from Hammer and Lyndersen,34 Perry and Weber,35 Small et al.,36 Dean et al.,37 Perry’s handbook,38 and Evans et al.39 They also performed their own experiments using two different methods, direct gas chromatography (DGC) and solid sorbent trapping/thermal desorption gas chromatography (STTD). It should be noted that the vapor pressures of TBP from Hammer and Lyndersen, Perry and Weber, Small et al., Dean, and Perry’s handbook were not directly measured using TBP. They utilized the DGC method to determine the retention times of TBP and dibutyl phthalate.33 The vapor pressures of TBP were then calculated using the retention times and previously reported vapor pressures of dibutyl phthalate at elevated temperature (170− 200 °C). Evans et al. measured the boiling points of TBP at three different pressures (50, 100, and 150 Torr) and fitted the data to another equation to obtain the vapor pressure of TBP at high temperatures.39 These experimental points were then extrapolated to the lower temperature range (0−40 °C) using linear regression. These kinds of analyses might lead to large errors due to few experimental points (three temperatures and three pressures) and the uncertainty of the extensive extrapolation (from 170−200 °C to 0−40 °C). However, the TBP vapor pressure determined using the STTD method by Skene et al.33 was measured directly using pure TBP at the lower range of temperatures usually found in the work place (0−80 °C), thus avoiding the uncertainties of extrapolation. The fitted value of the enthalpy of vaporization of TBP from this set is 22.67 kJ/mol. 2.3. EDM. The EDM the of TBP molecule, pD, was calculated by averaging over the instantaneous EDM of all molecules in the system pi =

⎧ ⎛ 1 1 ⎞⎫ I(g ) = I0 exp⎨−D(γδg )2 ⎜Δ − δ − τ ⎟⎬ ⎝ ⎩ 3 2 ⎠⎭

where I(g) and I0 represent integrated peak intensities in the presence and absence of gradient pulses of amplitude g. δ is the gradient duration, γ the gyromagnetic ratio of the nucleus observed, Δ the separation between gradient pulse pairs, and τ the time allowed for gradient recovery before the next pulse. In the TBP experiments, τ was set to 200 μs, and values of Δ and δ were manually optimized to 75 and 3.2 ms, respectively. On the basis of comparisons of repeated measurements as well as

∑ qa(ria − riCOM) a

pD =

⎛1 ⎜⎜ ⎝N

N

∑ i=1

1 d lim ⟨∑ |r(t ) − r(0)|2 ⟩ 6 t →∞ dt i = 1

⎞1/2

pi2 ⎟⎟ ⎠

C

DOI: 10.1021/jp510365c J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B error in the estimation of the center of intensity of the signals in indirect dimension, we estimate the overall accuracy of the diffusion coefficient reported here to be ±25%. 2.5. Self-Association of TBP in the n-Dodecane Mixture Study. 2.5.1. PMF Using the Weighted Histogram Analysis Method (WHAM). In order to study the selfassociation of TBP in n-dodecane, we first performed a PMF study using the umbrella sampling technique to establish the criteria for identifying TBP dimer and larger aggregates. The PMF study was carried out along two different reaction coordinates, distance and angle, as depicted in Figure 2: (1)

Table 1. Different Systems Studied volume fraction of TBP (%)

concentration (M)

# TBP molecules

# n-dodecane molecules

100 58 32.5 10 5 0

3.6 2.1 1.2 0.4 0.2 0

500 334 188 58 32 0

0 290 466 621 653 500

general self-association reactions, TBP dimers are formed by two TBP monomers coming together TBP + TBP ⇌ TBP2

The dimerization constant of TBP can be calculated by γTBP γTBP [TBP2] [TBP2] 2 K2 = × = × 2 2 2 γTBPγTBP [TBP][TBP] [TBP] γTBP where the square bracket denotes the concentration of the species and γ the activity coefficients on a molar scale. Similarly, a general reaction and equation can be derived for TBP aggregates TBPq − 1 + TBP ⇌ TBPq

Figure 2. Reaction coordinates for PMF study for two TBP molecules. (1) Distance in Å; (2) angle in degrees. Color scheme: red, oxygen atoms; gold, phosphorus atoms; cyan, carbon atoms; and white, hydrogen atoms.

Kq =

[TBPq] [TBP][TBPq − 1]

×

γTBP

q

γTBP γTBP q−1

It is common to assume that the activity coefficients in the organic phase remain constant and hence may be incorporated in the equilibrium constant, resulting in a conditional equilibrium constant, denoted Kq′. Using the criteria on distance and angle at which TBP molecules are energetically stable determined from the PMF study, we were able to calculate the concentrations of each TBP species present in the solutions, which in turn gave us the association constant of each species.

The distance (D) between the phosphorous (P) atom of one TBP molecule to the oxygen (O2) atom of the other TBP molecule ranging from 3.0 to 11.0 Å with a 0.5 Å increment. (2) The pseudo-dihedral angle (A) of the two PO2 double bonds of two separate TBP molecules rotating from 0 to 360° with 5° increment in the curved arrow direction. An angle of 0 or 360° represents the perfect antiparallel orientation, whereas 180° represents the perfect parallel orientation. Two TBP monomers were immersed in pure n-dodecane solvent and followed the simulation method described above for equilibration. During equilibration, position constraints were put on the P and O2 atoms so that the monomers were 10.5 Å apart and adapted the antiparallel orientation. During production runs, position restraints were placed on the P and O2 atoms to impose the reaction coordinates. The pseudodihedral angle was rotated from the initial 0 to 360° for each distance ranging from 3.0 to 11.0 Å. A harmonic force constant of 7 kcal/mol/Å2 was used for distance restraint, and 400 kcal/ mol/rad2 was used for angle restraint. A total of 1152 windows were simulated for 2 ns each. 1D and 2D PMFs were then calculated using the WHAM.43 2.5.2. RDFs and Self-Association Constants from MD Simulations of TBP/n-Dodecane Mixtures. In this study, we conducted long MD simulations of four TBP/n-dodecane mixtures at different concentrations ranging from 0.2 to 2.1 M, as shown in Table 1 at a constant temperature of 25 °C. At equilibrium, we calculated the densities of different mixtures and the RDFs between various atoms of different molecules by using a bin width of 0.25 Å and a cutoff distance of 32 Å (about half of the simulation box length), averaged over the last 20 ns of simulations. Also at equilibrium, we monitored the association between different TBP molecules. Following the

3. RESULTS AND DISCUSSION 3.1. Reparameterization of TBP and n-Dodecane Molecules. Similar to previous studies,18,19 which employed the OPLS-AA force field for n-dodecane molecules, we observed the same liquid-to-gel phase transition behavior from using the NPT ensemble and the default GAFF at a temperature of 298 K and pressure of 1 bar. Liquid-to-gel phase transition occurs when n-dodecane molecules initialized from a randomized liquid state (Figure 3A) interact with one another to form a much more ordered structure (Figure 3B). The density profile of the default GAFF showed pseudo-equilibrium in density at 2 and 3 ns, then increasing for the next 12 ns (Figure 3D). Thus, if a simulation is performed for only a short period of time, one might not observe this crystallization behavior of n-dodecane, leading to inaccurate results. Using our optimized parameter (OP) set without changing the default SFs, the liquid-to-gel phase transition of n-dodecane was completely eliminated, as showed by Figure 3C,D. The ndodecane density from this OP set remained around the experimental value of 0.74 g/cm3 over the whole simulation duration, and the molecules were still distributed randomly, representing the liquid phase (Figure 3C). D

DOI: 10.1021/jp510365c J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Table 3. Comparison of Densities, Enthalpies of Vaporization, and EDM of TBP and n-Dodecane Using the Default GAFF, the OPs Derived from This Studya, and Experimental Data ΔHvap (kJ/mol) TBP GAFF OP expt. n-Dodecane GAFF OP expt.

density (g/cm3)

25.63 ± 0.12 24.08 ± 0.04 22.67b

0.995 ± 0.012 0.971 ± 0.007 0.973c

17.14 ± 0.08 14.98 ± 0.04 14.70d

0.763 ± 0.0070f 0.744 ± 0.0005 0.745e

dipole (D) 2.88 ± 0.32 3.58 ± 0.34 3.51 ± 0.23c

a

Errors were estimated by a block average (6 blocks each 1 ns in length for the TBP molecule and 10 blocks each 1 ns in length for the n-dodecane molecule). bExperimental values are derived from Skene et al.33 cExperimental values are from ref 22. dExperimental values are from Haynes.44 eExperimental values are from Yaws.45 fThe starting density of the system was 0.742 g/cm3 and gradually increased to 0.763 g/cm3 at 10 ns with a positive slope.

Figure 3. Snapshots of a n-dodecane system starting as (A) an initial randomized configuration and ending as (B) an ordered configuration using GAFF or (C) a disordered configuration using the OPs. (D) The density profile of two simulations using GAFF and OPs is shown as a function of time.

The OP sets (Table 2) produced superior results when compared to experimental data as opposed to the default GAFF

density of pure TBP at four different temperatures ranging from 16.6 to 121.6 °C.46 Burgur recognized the linear relationship between the temperature and TBP’s density and reported an expansion coefficient of 0.000933 g/cm3/°C.22 He also recommended the value of 0.9727 ± 0.00004 g/cm3 for TBP density at room temperature after weighting the data at different temperatures from many previous studies.22 A comparison of simulated density of pure TBP at different temperatures with the calculated value is summarized in Table 4. We simulated pure TBP at five different temperatures from 20 to 100 °C using the NPT ensemble. The errors between the simulated and calculated values are less than 0.5%. 3.2. SDC. The MSD profiles of the system size of N = 500 molecules of the pure TBP and n-dodecane systems show two time regimes that correspond to different simulation ensembles, as shown in Figure 4. The first 15 ns of both systems were simulated under the NPT ensemble to equilibrate the systems’ densities, temperatures, and pressures, followed by 85 and 40 ns of the NVE ensemble for TBP and n-dodecane, respectively, to avoid introducing effects into the systems by employing the artificial thermostat and barostat. Indeed, the MSD slope under the NPT ensemble is remarkably different than that under the NVE ensemble. Such a difference between MSD slopes of the NPT part and the NVE part clearly demonstrate the effects of the thermostat and barostat. By relying on these regulators, the SDC might be underestimated as in the case of n-dodecane19 and overestimated as in the case of the TBP systems.15 It should be noted that when we performed the NVE simulations, special care was given to conserve the systems energy. In this case, we used an Ewald coefficient of 0.35 and a SHAKE tolerance of 10−7, which are more stringent than the default Amber values (0.275 for the Ewald coefficient and 10−5 for the SHAKE tolerance) to ensure that the system total energy drift is within the acceptable range for the GPU version of PMEMD.27 As a result, the temperatures deviate less than 2.5 K from 298 K over the 55 ns of simulation time for n-dodecane and over 100 ns of simulation time for TBP. This change in temperature should not have a significant impact on the reported SDC. The average pressure under the NVE ensemble is P = −6.9 ± 9.2 bar for TBP and P = 20.1 ± 11.4 bar for ndodecane. Even though these pressures fluctuated over a wider

Table 2. LJ Parameters and Partial Electric Charges for TBP and n-Dodecane Molecules Used in This Study atom type TBP O2 P O Ca Cb Cc Cd Ha Hb Hc Hd n-Dodecane Ca Cb Cc Ha Hb Hc

partial charge (q)

σ (Å)

ε (kcal/mol)

−0.8033 1.5955 −0.8033 0.1677 −0.0834 −0.0811 −0.0921 0.0454 0.0564 0.0417 0.0348

3.4945a 3.7418 3.0000 3.3997 3.3997 3.3997 3.3997 2.4714 2.6495 2.6495 2.6495

0.050a 0.2000 0.1700 0.1094 0.1094 0.1094 0.1094 0.0157 0.0157 0.0157 0.0157

−0.0921 −0.0804 −0.0794 0.0317 0.0387 0.0397

3.1324a 3.1324a 3.1324a 2.6495 2.6495 2.6495

0.1444a 0.1444a 0.1444a 0.0157 0.0157 0.0157

a

Values were modified from the original GAFF. Partial charges were obtained using antechamber package and AM1-BCC charge model.

for both molecules, as seen in Table 3. Table 3 shows the results on enthalpies of vaporization and densities of TBP and n-dodecane plus the EDM of TBP that were calculated using GAFF and our OPs. The differences between the OP and literature values are well within experimental errors. To our knowledge, we are the first group to report the value of the heat of vaporization of TBP via computational methods. Moreover, the densities of TBP and n-dodecane, the heat of vaporization of n-dodecane, and the dipole moment of TBP agreed well with previous optimization studies.15,18,19 The OP sets also produced a similar density of pure TBP as a function of temperatures. Vogel and Cowan measured the E

DOI: 10.1021/jp510365c J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B Table 4. Simulated Density of TBP at Different Temperatures temperature (°C) 20 25 50 75 100

simulated temperature (°C) 20.0 25.0 50.0 75.0 100

± ± ± ± ±

0.03 0.03 0.03 0.05 0.03

simulated pressure (bar) 0.98 0.98 0.98 0.96 0.98

± ± ± ± ±

0.09 0.08 0.19 0.15 0.11

simulated density (g/cm3) 0.976 0.971 0.948 0.924 0.901

± ± ± ± ±

0.0007 0.0007 0.0008 0.0004 0.0002

experimental valuea (g/cm3)

error (%)

± ± ± ± ±

0.15 0.22 0.27 0.19 0.25

0.977 0.973 0.950 0.926 0.903

0.0007 0.0004 0.0007 0.0007 0.0007

a Experimental values were obtained from the density of pure TBP at 25 °C of 0.9727 ± 0.00004 g/cm3 and the reported expansion coefficient of 0.000933 g/cm3/°C

Figure 4. MSD profiles of TBP and n-dodecane molecules (N = 500 molecules). Figure 6. (A) MSD of TBP as a function of time of four systems: blue, N = 200 molecules and L = 4.5 nm; red, N = 350 molecules and L = 5.4 nm; green, N = 500 molecules and L = 6.1 nm; purple, N = 1000 molecules and L = 7.7 nm. (B) SDC of TBP as a function of the inverse box length 1/L of the four systems.

range, they should not affect the SDC, as demonstrated by Jorgensen and co-workers, who studied the effect of pressure and temperature on the SDC of water and found that pressure plays such a minimal role on the SDC.47 The linear system size dependence of the SDC of ndodecane molecules, as observed previously by Yeh and Hummer,41 was confirmed, as shown in Figure 5. Figure 5A

independent of the system size, and significantly smaller than the experimental value of 2.29 × 10−6 cm2/s,21,22 and the value we obtained using 1H DOSY (2.8 ± 0.7 × 10−6 cm2/s; see Figures S3 and S4, Supporting Information). We suspect that it is due to the complexity of the TBP molecule that possesses both polar and nonpolar regions. We are in the process of parametrizing a polarizable force field in order to better represent this particular molecule. However, the uncorrected SDC of TBP, Ds = 2.1 × 10−7 cm2/s, is comparable to a reported value from a previous study.15 3.3. TBP/n-Dodecane Mixtures. 3.3.1. Mass Density. Different mixtures of TBP in n-dodecane with concentrations ranging from 0.2 to 2.1 M were simulated using the optimized force fields parametrized, as discussed earlier. The densities of different mixtures were calculated using the NPT ensemble and compared favorably, as showed by the small error in Table 5, with the results of Tian and Liu, who experimentally determined the densities of the binary mixtures over the entire composition range using a high-precision vibrating-tube digital density meter.49 3.3.2. RDFs. The RDFs between atoms belonging to different molecules were computed using a bin width of 0.25 Å and a cutoff distance of 32 Å (about half of the simulation box length), averaged over the last 20 ns of simulations. The RDFs can be used to gain insight into the structure of TBPs in the liquid phase as well as the organization of TBPs within the ndodecane solvent. In this study, we calculated the RDFs of O2−O2, O2−P, and P−P atom pairs at 0.2 M TBP concentration in n-dodecane solution, as shown in Figure 7. The O2−O2 RDF shows a distinctive peak at 4.9 Å (Figure 7A); the O2−P RDF shows a sharper peak at 5.1 Å (Figure 7B); and the P−P RDF possesses the narrowest peak at 5.9 Å (Figure 7C). Together with the PO double bond length of 1.48 Å, the O2−O2 distance of 4.9 Å (Figure 7A) and the O2−P distance of 5.1 Å (Figure 7B) form a right angle, suggesting the dipole−dipole orientation of

Figure 5. (A) MSD of n-dodecane as a function of time of three systems: blue, N = 200 molecules and L = 4.2 nm; red, N = 350 molecules and L = 5.1 nm; green, N = 500 molecules and L = 5.8 nm. (B) SDC of n-dodecane as a function of the inverse box length 1/L of the three systems.

shows that the MSDs grow linearly with time after the equilibration period of 15 ns at NPT and more steeply with increasing system size. Figure 5B shows the apparent SDC of the n-dodecane molecule as a function of the inverse box length (1/L). By extrapolating the linear fit to the apparent SDC profile, a corrected SDC of Dscorr = 7.24 × 10−6 cm2/s was obtained for an infinite system of n-dodecane molecules at room temperature and ambient pressure. This new corrected SDC of n-dodecane molecule, even though still deviates 17% from the experimental value of 8.71 × 10−6 cm2/s,48 is the best estimation compared to other studies by the Khomami group and Bockmann group, who reported values of 4.64 ± 0.234 × 10−6 cm2/s using SF = 0.4 and 5.29 ± 0.13 × 10−6 cm2/s, respectively. However, the same correction study for system size dependence did not work when applied to TBP molecules (Figure 6). We found that the SDC of TBP was constant, F

DOI: 10.1021/jp510365c J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B Table 5. Simulated Density Values of TBP/n-Dodecane Mixtures at 25°C TBP concentration (M) 2.1 1.2 0.4 0.2

simulated temperature (°C) 25.00 25.01 24.98 25.00

± ± ± ±

0.01 0.03 0.01 0.04

simulated pressure (bar) 0.50 0.28 −0.50 0.78

± ± ± ±

1.74 1.57 0.78 1.65

simulated density (g/cm3)

experimental value (g/cm3)

error (%)

0.870 0.814 0.766 0.758

0.871 0.815 0.770 0.761

0.21 0.20 0.49 0.66

M−1 for TBP in n-dodecane without explicitly specifying the concentration range.52 They used the same mathematical model for both solvents with the assumption that there were only TBP monomers and dimers present in the solution at all studied concentrations. In our self-association study, following the PO dipole− dipole interaction suggestion, we performed umbrella sampling simulations in which two TBP molecules were arranged in a fixed antiparallel orientation and the distance between the P atom of one molecule and the O atom of the other one was varied. Figure 8A shows that it is energetically unfavorable for

Figure 7. RDFs of (A) O2−O2 atoms, (B) O2−P atoms, and (C) P− P atoms on any two different TBP molecules at 0.2 M TBP concentration in n-dodecane solution. (D) Proposed models of a TBP dimer in the antiparallel orientation; the butyl tails are omitted for clarity.

two TBP molecules where the PO double bonds are positioned antiparallel to each other, as shown in Figure 7D. This proposed orientation is in agreement with the results from our PMF study in which the global free-energy minimum exists for the TBP dimers with the antiparallel orientation and are ∼5.5 Å apart, as discussed below. This antiparallel orientation is owed to the presence of the large butyl groups, thus making it difficult for two TBP molecules to approach each other from the hydrophobic sides, and the fact that the polar heads PO prefer to interact with one another in order to minimize the system energy. The large butyl groups also contribute to the existence of only one coordination shell (CN) for TBP molecules at a low concentration of 0.2 M. The first-shell CN of TBP was computed by integrating the RDFs of the P−P atom pair up to the minimum of the first peaks. The CN of TBP is 0.86 at 0.2 M. The RDFs of all atom pairs approach the g(r) value of 1 at around 13 Å. This correlation length of TBP molecules implies that they stop having effects on theirs neighbors when they are 13 Å apart. This is the reason that we picked the longest distance for our PMF study to be 11 Å. 3.3.3. Self-Association of TBP in n-Dodecane Mixtures. It has been suggested that the PO dipole−dipole interaction is responsible for the self-association of TBP molecules based on previous studies.50−53 Different techniques such as vapor pressure osmometry53 and IR spectroscopy51,52 had been utilized for the determination of the TBP dimerization constant as a function of its environments. Petkovic reported a value of 2.9 ± 0.1 M−1 for the dimerization of TBP in n-hexane with concentrations from 0.2 to 1.0 M51 and a value of 2.6 ± 0.1

Figure 8. 1D PMF of a TBP dimer as a function of (A) distance while the pseudo-dihedral angle is fixed at 360° and (B) the pseudo-dihedral angle while the distance is fixed at 5.25 Å.

two TBP molecules, even with the antiparallel orientation, to be too close to each other around 3 Å. A global free-energy minimum is observed at 5.25 Å for TBP dimers, in agreement with the peaks in the RDF study of the previous section. Next, we fixed the distance between the two molecules to be at the optimal distance of 5.25 Å and let one of the TBP molecules rotate a full period (0−360°) to study the effect of orientation on the stability of TBP dimers. TBP dimers clearly favor the antiparallel arrangement, as can be seen by the low free energy of the pseudo-dihedral angles in the range of 0−60 or 300− 360° (Figure 8B). However, these results did not provide a complete picture of the free-energy map for TBP dimers over the wide range of distances and angles. Therefore, a 2D PMF study was carried out in which we sampled 1152 windows of different distances and angles for 2 ns each. The contour map of the free energy along reaction coordinates of the pseudo-dihedral angles and distances between two TBP molecules is presented in Figure 9A. There is the expected symmetry in the angle reaction coordinate as shown by the blue regions. When two TBP molecules are close to one another (D < 4.0 Å), it is difficult to G

DOI: 10.1021/jp510365c J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 10. Self-association of TBP in 0.2 M TBP/n-dodecane solution profiles over 180 ns of simulation showing the percentage of different TBP species: monomer, dimer, and trimer.

ns, there is an average of 65% of TBP molecules remaining as monomers with a standard deviation of 6.4%. The large fluctuation in the number of monomers present in the system is indicative of the sampling process where TBP molecules come together and break apart in order to minimize the system’s energy. During the last 40 ns of the self-association profile, the average amount of TBP monomers reduced to 64% with a lower standard deviation of 4.6%, demonstrating that the system had reached equilibrium. The average number of TBP trimers in the mixture is less than 1 (9%) at all times. With these values, the dimerization and trimerization constants were calculated, using the equations described in a previous section, to be 1.2 ± 0.03 and 2.2 ± 0.2 M−1, respectively. In order to compare our simulated dimerization constant with experimental measurement,52 we made the same assumption as the experimental study that there were only monomers and dimers existing in the solution. The dimerization constant then was determined to be 2.2 ± 0.5 M−1. This value matches well with the experimental value of 2.6 ± 0.1 M−1 reported by Petkovic, who performed the IR spectroscopic technique.52

Figure 9. (A) The 2D PMF profile along the reaction coordinate distance and angle obtained from the umbrella sampling method and WHAM analysis of two TBP molecules. (B) Two TBP molecules adopt the antiparallel orientation and form a stable dimer with low free energy. (C) Two TBP molecules assume parallel arrangement and possess high free energy due to steric effect. The arrows represent the direction of each TBP’s permanent dipole moment.

sample structures in the regions of angles around 100−300°, as demonstrated by the red regions that correspond to high freeenergy values. They can only assume the antiparallel orientation due to steric effect, caused by the bulky butyl tails, and possibly the high electrostatic interaction energy of like-sign charges. As they move farther away, their free energies minimize but with a clear preference for the antiparallel (Figure 9B) over the parallel orientation (Figure 9C). The free-energy minima are observed for distances above 4.0 Å with angles in the range of 0−60 and 300−360°. Specifically at distances around 5.0−6.0 Å, the global free-energy minimum is observed over the angle range of 0−30 and 330−360°. Recently, Khomami and co-workers also studied the physical properties and molecular structure of TBP/n-dodecane mixtures17 using various AMBER and OPLS force fields for TBP and OPLS-AA for n-dodecane with different SFs for relatively short MD simulations. They computed the angular distribution function over a range of fixed distances between 4.0 and 7.0 Å and found that this angular distribution peaks at around 40° less than the perfect alignment in which the two dipoles of a TBP dimer are in antiparallel orientation using the OPLS-DFT/RHF model for TBP. Moreover, when using the OPLS-MNDO model for TBP, they also found that the angular distribution peaks at around 90°, which belongs to a perpendicular orientation. These results disagree with our RDF and PMF results. We suspect that one of the reasons for this disagreement is related to the length of MD simulations that they conducted in comparison to our relatively long simulations and the exhaustively conformational search techniques of umbrella sampling and WHAM that we used to calculate PMFs. The distances and angles of the free-energy minima were used as criteria for the determination and quantification of TBP aggregates in a 0.2 M TBP/n-dodecane concentration, as shown in Figure 10. In the beginning of the simulation, it took about 10 ns for the system to reorganize and move away from the initial coordination generated by Packmol. For the next 130

4. CONCLUSIONS In this work, we have successfully parametrized new parameter sets for TBP and n-dodecane molecules using the AMBER force field. Parameterization was performed on the LJ potential parameters to match the experimental heat of vaporization and densities that are dependent upon the temperature. The dipole moment and SDC were used as means to validate the OP sets. The values resulting from our OPs show better accuracy than the ones from the default GAFF when compared with experimental values. More importantly, the liquid-to-gel phase transition of n-dodecane observed with the default GAFF (or OPLS-AA18,19) was completely eliminated. This demonstrates that instead of combining multiple force fields for different molecules within the same system,17 one can perform rigorous, self-consistent, and systematic parametrization by refining the LJ potential parameters on just one force field such as GAFF. The encouraging results of our parametrization study for two very dissimilar molecules such as TBP and n-dodecane suggests that this parametrization method is not limited to only this study but can be applied to any molecules of interest. Within this method, the NPT ensemble should be employed to correctly equilibrate the system’s density, and a long equilibration period should be adopted to avoid pseudoequilibrium states as observed with the n-dodecane density in this study. We also examined the structure of TBPs in a TBP/ndodecane mixture by determining the RDFs at equilibrium of MD simulations on a mixture and calculating PMFs using H

DOI: 10.1021/jp510365c J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

U.S. Department of Energy for funding the work through the Nuclear Energy University Program, NEUP Contract No. 120569. Q.N.V. acknowledges support from a Graduate Research Fellowship from the National Science Foundation (DGE-1321846). The U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences, Division of Chemical Sciences, Geosciences and Biosciences funded the work performed by L.X.D.

WHAM. The results from RDF study indicate that the most probable orientation of two TBP molecules in the TBP/ndodecane mixtures is antiparallel. This observation agrees with the PMF study, which showed a clear preference of the antiparallel orientation over the parallel one. The criteria extracted from the 2D PMF free-energy minima were used for the calculation of the dimerization constant of TBP in ndodecane. To our knowledge, we are the first group that is able to obtain a dimerization constant of 2.2 ± 0.5 M−1 from MD simulations, and this value compares favorably with the reported value of 2.6 ± 0.1 M −1 determined by IR spectroscopy.52 For the SDC of the n-dodecane molecule, we observed that it was necessary to perform the correction to remove the system size dependency, as recommended by the work of Yeh and Hummer.41 However, this correction was not applicable for the SDC of the TBP molecule due to its complexity, where the TBP molecule possesses both polar and nonpolar regions. We are in the process of parametrizing a polarizable force field for this particular molecule in order to better present TBP in the future. We also recommend utilizing the microcanonical ensemble (NVE) when calculating the SDC. SDC is a dynamic property, which would be perturbed by the use of a thermostat (in NPT or NVT ensembles) because temperature is controlled by the scaling of the molecules’ velocities. The effect of the thermostat on the SDC value was clearly demonstrated throughout our study. Additionally, extended simulation durations are required to achieve the linear Einstein relation for the MSD over time. Regarding the self-association of TBP in TBP/n-dodecane mixtures, we performed our study for the whole range of concentrations described in the beginning of this work but only reported the result at the low concentration of 0.2 M. To our knowledge, there are no experimental studies reported for the self-association of TBP in n-dodecane solutions at high concentrations due to the complication of forming larger aggregate sizes. We plan to conduct our own IR experiments that will include models for higher oligomer sizes of TBP. The results from this future work will greatly aid the understanding of the solvent extraction process at high TBP concentrations.

■

ABBREVIATIONS

■

REFERENCES

SX, solvent extraction; PUREX, Plutonium Uranium Reduction Extraction; TBP, tri-n-butyl phosphate; MD, molecular dynamics; SDC, self-diffusion coefficient; LJ, Lennard-Jones; GAFF, general AMBER force field; RDF, radial distribution function; PMF, potential mean force; NPT, isothermal− isobaric; NVE, microcanonical; DGC, direct gas chromatography; STTD, sorbent trapping/thermal desorption; EDM, electric dipole moment; MSD, mean-squared displacement; WHAM, weighted histogram analysis method; OP, optimized parameter; CN, coordination number

(1) Problems Concerning the Accumulation of Isolated Plutonium; Bay, H., Ed.; IAEA-TECDOC-765; International Atomic Energy Agency: Vienna, Austria, 1994. (2) Sood, D. D.; Patil, S. K. Chemistry of Nuclear Fuel Reprocessing: Current Status. J. Radioanal. Nucl. Chem. 1996, 203, 547−573. (3) Lanham, W. B.; Runion, T. C. Purex Process for Plutonium and Uranium Recovery; ORLN-497; DOE Report; Department of Energy: Washington, DC, 1949. (4) Schulz, W. W.; Burger, L. L.; Navratil, J. D.; Bender, K. P. Science and Technology of Tributyl Phosphate; CRC Press: Boca Raton, FL, 1990; Vol. 3. (5) Zilberman, B. Y.; Fedorov, Y. S.; Kopyrin, A. A.; Arkhipov, S. A.; Blazheva, I. V.; Glekov, R. G. Extraction of U(IV) and U(VI) with 30% Tributyl Phosphate Under Conditions of Formation of the Second Organic Phase. Radiochemistry 2001, 43, 172−176. (6) Chiarizia, R.; Nash, K. L.; Jensen, M. P.; Thiyagarajan, P.; Littrell, K. C. Application of the Baxter Model for Hard Spheres with Surface Adhesion to SANS Data for the U(VI)−HNO3, TBP−n-Dodecane System. Langmuir 2003, 19, 9592−9599. (7) Chiarizia, R.; Jensen, M. P.; Borkowski, M.; Ferraro, J. R.; Thiyagarajan, P.; Littrell, K. C. Third Phase Formation Revisited: the U(VI), HNO3−TBP, n-Dodecane System. Solvent Extr. Ion Exch. 2003, 21, 1−27. (8) Chiarizia, R.; Jensen, M. P.; Rickert, P. G.; Kolarik, Z.; Borkowski, M.; Thiyagarajan, P. Extraction of Zirconium Nitrate by TBP in nOctane: Influence of Cation Type on Third Phase Formation According to the “Sticky Spheres” Model. Langmuir 2004, 20, 10798−10808. (9) Plaue, J.; Gelis, A.; Czerwinski, K.; Thiyagarajan, P.; Chiarizia, R. Small-Angle Neutron Scattering Study of Plutonium Third Phase Formation in 30% TBP/HNO3/Alkane Diluent Systems. Solvent Extr. Ion Exch. 2006, 24, 283−298. (10) Beudaert, P.; Lamare, V.; Dozol, J. F.; Troxler, L.; Wipff, G. Theoretical Studies on Tri-n-Butyl Phosphate: Md Simulations in Vacuo, in Water, in Chloroform, and at a Water/Chloroform Interface. Solvent Extr. Ion Exch. 1998, 16, 597−618. (11) Burgard, M.; Wipff, G. TBP at the Water−Oil Interface: The Effect of TBP Concentration and Water Acidity Investigated by Molecular Dynamics Simulations. J. Phys. Chem. B 2001, 105, 11131− 11141. (12) Ye, X.; Cui, S.; de Almeida, V.; Khomami, B. Interfacial Complex Formation in Uranyl Extraction by Tributyl Phosphate in Dodecane Diluent: A Molecular Dynamics Study. J. Phys. Chem. B 2009, 113, 9852−9862.

ASSOCIATED CONTENT

S Supporting Information *

Force field reparameterization of TBP and n-dodecane and DOSY spectra from NMR experiements. This material is available free of charge via the Internet at http://pubs.acs.org.

■

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Tel.: 949-824-6589. Author Contributions

The manuscript was written through contributions of all authors. All authors have given approval to the final version of the manuscript. All authors contributed equally. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We would like to give special thanks to Phil Dennison (UCI NMR Facility Director) for helping with the NMR experiments and the High Performance Computing (HPC) Cluster at UCI for computational resources. The authors wish to thank the I

DOI: 10.1021/jp510365c J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B (13) Ye, X.; Smith, R. B.; Cui, S.; de Almeida, V.; Khomami, B. Influence of Nitric Acid on Uranyl Nitrate Association in Aqueous Solutions: a Molecular Dynamics Simulation Study. Solvent Extr. Ion Exch. 2010, 28, 1−18. (14) Ye, X.; Cui, S.; de Almeida, V. F.; Hay, B. P.; Khomami, B. Uranyl Nitrate Complex Extraction Into TBP/Dodecane Organic Solutions: A Molecular Dynamics Study. Phys. Chem. Chem. Phys. 2010, 12, 15406−15409. (15) Hay, B. P.; Cui, S.; Khomami, B.; de Almeida, V. F.; Ye, X. Molecular Dynamics Simulation of Tri-n-butyl-phosphate Liquid: A Force Field Comparative Study. J. Phys. Chem. B 2012, 116, 305−313. (16) Ye, X.; Cui, S.; de Almeida, V. F.; Khomami, B. Molecular Simulation of Water Extraction Into a Tri-n-butylphosphate/ nDodecane Solution. J. Phys. Chem. B 2013, 117, 14835−14841. (17) Cui, S.; de Almeida, V. F.; Khomami, B. Molecular Dynamics Simulations of Tri-n-butyl-phosphate/n-Dodecane Mixture: Thermophysical Properties and Molecular Structure. J. Phys. Chem. B 2014, 118, 10750−10760. (18) Siu, S. W. I.; Pluhackova, K.; Böckmann, R. A. Optimization of the OPLS-AA Force Field for Long Hydrocarbons. J. Chem. Theory Comput. 2012, 8, 1459−1470. (19) Ye, X.; Cui, S.; de Almeida, V. F.; Khomami, B. Effect of Varying the 1−4 Intramolecular Scaling Factor in Atomistic Simulations of Long-Chain n-Alkanes with the OPLS-AA Model. J. Mol. Model. 2012, 19, 1251−1258. (20) Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A. Development and Testing of a General Amber Force Field. J. Comput. Chem. 2004, 25, 1157−1174. (21) Johnson, W.; Dillon, R. Physical Properties of TributylphosphateDiluent Solutions; HW-29086; 1953. (22) Schulz, W. W.; Navratil, J. D.; Talbot, A. E. Science and Technology of Tributyl Phosphate; CRC Press: Boca Raton, FL, 1984; Vol. 1. (23) Jakalian, A.; Bush, B. L.; Jack, D. B.; Bayly, C. I. Fast, Efficient Generation of High-Quality Atomic Charges. AM1-BCC Model: I. Method. J. Comput. Chem. 2000, 21, 132−146. (24) Jakalian, A.; Jack, D. B.; Bayly, C. I. Fast, Efficient Generation of High-Quality Atomic Charges. AM1-BCC Model: II. Parameterization and Validation. J. Comput. Chem. 2002, 23, 1623−1641. (25) Wang, J.; Wang, W.; Kollman, P. A. Automatic Atom Type and Bond Type Perception in Molecular Mechanical Calculations. J. Mol. Graphics Modell. 2006, 25, 247−260. (26) Darden, T. A.; T E Cheatham, I.; Simmerling, C. L.; Wang, J.; Duke, R. E.; Luo, R.; Walker, R. C.; Zhang, W.; Merz, K. M.; Roberts, B.; et al. Amber 12; University of California, San Francisco, CA, 2012. (27) Götz, A. W.; Williamson, M. J.; Xu, D.; Poole, D.; Le Grand, S.; Walker, R. C. Routine Microsecond Molecular Dynamics Simulations with AMBER on GPUs. 1. Generalized Born. J. Chem. Theory Comput. 2012, 8, 1542−1555. (28) Le Grand, S.; Götz, A. W.; Walker, R. C. SPFP: Speed Without Compromise  A Mixed Precision Model for GPU Accelerated Molecular Dynamics Simulations. Comput. Phys. Commun. 2013, 184, 374−380. (29) 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− 2164. (30) Qiao, B.; Demars, T.; Olvera de la Cruz, M.; Ellis, R. J. How Hydrogen Bonds Affect the Growth of Reverse Micelles Around Coordinating Metal Ions. J. Phys. Chem. Lett. 2014, 5, 1440−1444. (31) Dang, L. X.; Chang, T.-M. Molecular Dynamics Study of Water Clusters, Liquid, and Liquid−Vapor Interface of Water with ManyBody Potentials. J. Chem. Phys. 1997, 106, 8149. (32) Dang, L. X. Intermolecular Interactions of Liquid Dichloromethane and Equilibrium Properties of Liquid−Vapor and Liquid− Liquid Interfaces: A Molecular Dynamics Study. J. Chem. Phys. 1999, 110, 10113.

(33) Apelblat, A.; Hornik, A. Application of Vapor Pressure Osmometry to the Study of Uranyl Nitrate and Thorium Nitrate Complexes with Tributyl Phosphate. Isr. J. Chem. 1969, 7, 45−48. (34) Hammer, E.; Lydersen, A. L. The Vapour Pressure of Di-nbutylphthalate, Di-n-butylsebacate, Lauric Acid and Myristic Acid. Chem. Eng. Sci. 1957, 7, 66−72. (35) Perry, E. S.; Weber, W. H. Vapor Pressures of Phlegmatic Liquids. II. High Molecular Weight Esters and Silicone Oils. J. Am. Chem. Soc. 1949, 71, 3726−3730. (36) Small, P. A.; Small, K. W.; Cowley, P. The Vapour Pressures of Some High Boiling Esters. Trans. Faraday Soc. 1948, 44, 810−816. (37) Dean, J. A.; Lange, N. A. Lange’s Handbook of Chemistry, 13th Ed.; McGraw-Hill: New York, 1984. (38) Perry’s Chemical Engineers’ Handbook, 6th Ed.; McGraw-Hill: New York, 1984. (39) Evans, D. P.; Davies, W. C.; Jones, W. J. CLXVII.  The Lower Trialkyl Orthophosphates. Part I. J. Chem. Soc. 1930, 1310. (40) Gereben, O.; Pusztai, L. System Size and Trajectory Length Dependence of the Static Structure Factor and the Diffusion Coefficient as Calculated from Molecular Dynamics Simulations: The Case of SPC/E Water. J. Mol. Liq. 2011, 161, 36−40. (41) Yeh, I.-C.; Hummer, G. System-Size Dependence of Diffusion Coefficients and Viscosities from Molecular Dynamics Simulations with Periodic Boundary Conditions. J. Phys. Chem. B 2004, 108, 15873−15879. (42) Wu, D. H.; Chen, A. D.; Johnson, C. S. An Improved DiffusionOrdered Spectroscopy Experiment Incorporating Bipolar-Gradient Pulses. J. Magn. Reson., Ser. A 1995, 115, 260−264. (43) Grossfield, A. WHAM: The Weighted Histogram Analysis Method, version 2.0.9. http://membrane.urmc.rochester.edu/content/wham (2012). (44) Haynes, W. M. CRC Handbook of Chemistry and Physics, 95th ed.; CRC Press: Boca Raton, FL, 2014. (45) Yaws, C. L. In Yaws’ Critical Property Data for Chemical Engineers and Chemists; Yaws, C. L., Ed.; Knovel: New York, 2012. (46) Vogel, A. I.; Cowan, D. M. Physical Properties and Chemical Constitution. Part VII. Alkyl Sulphides, Disulphides, Sulphites, Sulphates, and Orthophosphates. J. Chem. Soc. 1943, 16. (47) Mahoney, M. W.; Jorgensen, W. L. Diffusion Constant of the TIP5P Model of Liquid Water. J. Chem. Phys. 2001, 114, 363. (48) Tofts, P. S.; Lloyd, D.; Clark, C. A.; Barker, G. J.; Parker, G. J.; McConville, P.; Baldock, C.; Pope, J. M. Test Liquids for Quantitative MRI Measurements of Self-Diffusion Coefficient in Vivo. Magn. Reson. Med. 2000, 43, 368−374. (49) Tian, Q.; Liu, H. Densities and Viscosities of Binary Mixtures of Tributyl Phosphate with Hexane and Dodecane from (298.15 to 328.15) K. J. Chem. Eng. Data 2007, 52, 892−897. (50) Miller, C. D.; Miller, R. C.; William Rogers, J. Phosphine Oxides. V. Intra- and Intermolecular Association. J. Am. Chem. Soc. 1958, 80, 1562−1565. (51) Petković, D. M. Dimerization of Tri-n-butyl Phosphate. Proceedings of the International Conference on Solvent Extraction Chemistry; Gothenburg: North Holland, Amsterdam, 1967; p 305. (52) Petković, D. M. Some Correlations of Trialkyl Phosphates Dimerization Constants. J. Inorg. Nucl. Chem. 1968, 30, 603−609. (53) Huang, C.-H.; Bautista, R. G. The Aggregation and Interactions of Tributyl Phosphate and Tricaprylmethylammonium Nitrate in Hexane by Osmometry. Sep. Sci. Technol. 2006, 19, 515−529.

J

DOI: 10.1021/jp510365c J. Phys. Chem. B XXXX, XXX, XXX−XXX