CuCl Complexation in the Vapor Phase : Insights from Ab Initio Molecular Dynamics Simulations

We investigated the hydration of the CuCl complex in HCl-bearing water vapor at 350C and a vapor-like fluid density between 0.02 and 0.09 g/cm using ab initio molecular dynamics (MD) simulations. The simulations reveal that one water molecule is strongly bonded to Cu(I) (first coordination shell), forming a linear [H2O-Cu-Cl] 0 moiety. The second hydration shell is highly dynamic in nature, and individual configurations have short life-spans in such low-density vapors, resulting in large fluctuations in instantaneous hydration numbers over a timescale of picoseconds.The average hydration number in the second shell (m) increased from ∼0.5 to ∼3.5 and the calculated number of hydrogen bonds per water molecule increased from 0.09 to 0.25 when fluid density (which is correlated to water activity) increased from 0.02 to 0.09 g/cm (fH2O 1.72 to 2.05). These changes of hydration number are qualitatively consistent with previous solubility studies under similar conditions, although the absolute hydration numbers from MDwere much lower than the values inferred by correlating experimental Cu fugacity with water fugacity.This could be due to the uncertainties in theMD simulations and uncertainty in the estimation of the fugacity coefficients for these highly nonideal “vapors” in the experiments. Our study provides the first theoretical confirmation that beyond-first-shell hydrated metal complexes play an important role in metal transport in low-density hydrothermal fluids, even if it is highly disordered and dynamic in nature.


Introduction
Understanding the speciation of metals in hydrothermal fluids is important as it controls metal mobility and mineral solubility in geological and man-made hydrothermal systems, such as those forming ore deposits and geothermal systems and those used in hydrometallurgy [1].Recent experimental and theoretical studies have shown that metals such as Cu and Mo can be transported in both liquid and vapor phases [2][3][4].Copper is transported as chloride species in many hydrothermal fluids [5].In high-density, reducedsulfur-poor brines, the dominant copper species is CuCl 2 [6][7][8][9][10], and the activity of the chloride ion is the main control on the solubility of Cu minerals at constant temperature () and pressure ().In contrast, in low-density vapors, there is increasing experimental evidence showing that, aside from the fugacity of HCl 0 , the solubility of Cu minerals is controlled primarily by fluid density, which correlates with the fugacity of water.Specifically, mineral solubility increases exponentially with increasing fluid density in HCl-bearing water vapor [3,11,12].
Enabled by advances in computing power and algorithms, molecular dynamics simulations have been applied increasingly to the study of metal speciation in hydrothermal liquids (e.g., Cu(I)-Cl -, [13,14]; Au(I)-HS -, [15]; Ag(I)-Cl -, [16,17]; Au(I)-HS -/OH -/S 3 -, [18,19]; Cu(I)-HS -/Cl -, [20]; Au(I)-Cl -, [14]; Pd(II)-Cl -/HS -, [21]; Zn(II)-Cl -/HS -, [22,23]; Ag(I)-HS -/OH -, [24]; Cu(I)-Ac -, [25]).These studies have demonstrated the importance of MD in understanding metal complexation at the molecular-level and, in some cases, have provided complementary information that helped resolve apparent discrepancies in the interpretation of experimental results [13,22].Molecular dynamics simulations have also been used to study metal transport in low-density aqueous fluids, namely, to determine the hydration behavior of Cu(I)-Cl and Au(I)-Cl species in fluids with a density as low as 0.1 g/cm 3 [14], and the speciation and hydration of Au(I)-Cl complex in CO 2 -rich vapors [26].Concurrently, Migdisov et al. [3] measured the solubility of Cu(I) chloride in HCl-bearing water vapor between 350 and 550 ∘ C, with a fluid density down to 0.013 g/cm 3 at 350 ∘ C.They found that the solubility of CuCl(s) increased exponentially as a function of water fugacity and proposed that this increase reflects an increasing number of hydration waters around the CuCl 0 species.They reported formation constants for hydrated CuCl⋅nH 2 O 0 clusters containing one to fourteen water molecules (n = [1][2][3][4][5][6][7][8][9][10][11][12][13][14]. The aim of this study is to conduct MD simulations under conditions similar to those of the experiments of Migdisov et al. [3], especially at conditions with significant change of Cu solubility (f H 2 O of 1.72-2.05,or fluid density 0.02-0.09g/cm 3 ; 350 ∘ C) to (1) verify the validity of the assumption that the increasing solubility of Cu with increasing water fugacity reflects a rapid change in hydration number; (2) to determine the nature and dynamics of the hydration sphere; and (3) to characterize the coordination structure of the dominant Cu species in HCl-bearing water vapor.

Method and Simulation Details
2.1.Ab initio Molecular Dynamics Simulations.Ab initio molecular dynamics simulations of CuCl complexes in the vapor phase were performed using the Car-Parrinello MD code "CPMD" version 3.17.1 [27].Details of the calculation procedures are given in Mei et al. [14]; the key elements are repeated here.Vanderbilt ultrasoft pseudopotentials in the CPMD package were used with plane-wave cutoffs of 25 Ry (340.14 eV), a time-step of 4 a.u.(0.097 fs), and a fictitious electron mass of 400 a.u.(3.644 × 10 -28 kg).The PBE exchange correlation-functional was used with a cut-off gradient correction of 5 × 10 −5 .Temperature was controlled by the Nosé thermostat for both the ions and electrons.All MD simulations were performed in the NVT ensemble at 350 ∘ C and run for >50 ps.The density of water was varied from 0.02 to 0.09 g /cm 3 in the simulations and the pressure () taken from the equations-of-state in the NIST database (Lemmon et al., 2000).The fugacity of water was calculated with HCh package [28] using water model of Kestin et al. [29].
The composition of the simulation box (1 Cu, 1 Cl, and 55 H 2 O) was chosen to represent the CuCl 0 complex in a vapor with low Cu and Cl concentrations.The formal concentrations of Cu and Cl are 1 molal, but since there were only one Cu and one Cl in the simulation box, there were no Cl-Cl or Cu-Cu interactions, which is consistent with the high dilution in vapors.In addition, the CuCl 0 species (Figure 1) was stable over the full simulation time in all the simulations.Hence, the selected model box reflects the results of experimental studies well, as they [3,11] have shown that a hydrated CuCl 0 complex is the dominant species under the simulated conditions.The simulation details are listed in Table 1.

Hydration Number and H-Bonds.
To obtain the time average for the ion-pairing and the hydration number, the radial distribution functions (RDF) and their integrals   [31], a Hbond between the donor (D) and acceptor O or Cl atom (A) would be present, if the D-H distance is less than 3.5 Å and the D-H-A angle is less than 70 ∘ .We also quantified the average number of H-bonds among bulk water molecules in the simulation box using the same geometrical approach.

Overall Geometry of the First-Shell CuCl Complex.
All MD simulations produced a linear first-shell complex, [CuCl(H 2 O)] 0 (Figure 1).The CuCl distance (2.05-2.07Å), Cu-O distance (1.91-1.93Å), and O-CuCl angle (162-166 ∘ ) at 350 ∘ C did not differ significantly over the density range of 0.02-0.09g/cm 3 (Table 1).The bond distance and angles are consistent with results of the in situ X-ray absorption spectroscopy (XAS) studies of Cu(I)-Cl and Cu(I)-HS complexes [6,12,32].The CuCl distance determined in these simulations is also similar to that determined by Mei et al. [14] for the linear CuCl 2 − complex for high-density (0.67-1.03 g/cm 3 ) solutions and higher ionic strength (4 m Cl − ) at 300 ∘ C.   2. The number of hydration water molecules based on the integral of the RDF increases with density from 2.7 (1.7, excluding the first-shell oxygen) at a density of 0.02 g/cm 3 , to 3.8 at 0.03 g/cm 3 , and to 4 at 0.08 and 0.09 g/cm 3 .The evolution of the instantaneous hydration number at different values of fluid density is shown in Figure 3.The dynamic hydration number displays a large variation as a function of simulation time, indicating that the stability of the system is lower than is the case in highdensity solutions [14].In other words, the second-shell water molecules of CuCl 0 in the vapor have short residence times; individual configurations have lifetimes ≤ 2 ps, similar, for example, to the lifetimes of sodium chloride ion pairs at 380 ∘ C, 0.55 g/cm 3 [33].

Water-Water Interaction.
Here we compare the environment of solvent water molecules (bulk water) to that of the   the CN between the 0.02-0.025g/cm 3 density group and the 0.03-0.09g/cm 3 density groups reflects the decrease in hydration with decreasing density.

Chloride Hydration. The distribution of solvent water
O around the Cl bound to Cu(I) is related to peaks at 3.5-3.6Å (Figure 6).The intensity of these RDF peaks is smaller than that for the O bound to solvent water O (Figure 5(b)), indicating weaker interaction between Cl and solvent water.This suggests that the water cluster surrounding the bound water contributes more to Cu hydration than the water cluster surrounding chloride.If we consider water clusters as acting as ligands in a manner similar to firstshell Cl -and H 2 O, the CuCl 0 complex can be written as where  is the number of water molecules surrounding Cl and  is the number of water molecules surrounding the bound water in the second shell.

Hydrogen
Bonds.An explanation for the strong interaction between solvent water and the [CuCl(H 2 O)] 0 complex is that water molecules form clusters that attach to metal complexes with the aid of strong hydrogen bonds (H-bonds) [34].The number of H-bonds per water molecule in the simulation boxes estimated from the bond geometry is listed in Table 2. Figure 7 shows four snapshots of the simulation for a density of 0.08 g/cm 3 .There H atoms play the roles of donor (D) and Cl or O plays the role of acceptor (A) in the H-bonds.Similar hydrogen bonding structure and the hydration of Cl -in hydrothermal NaCl solutions were also quantitatively addressed by experimental (IR and Raman) and computational (classical MD) approaches [35].They found that the feature of the NaCl structure has a considerable    as observed by previous studies for H-bonds (the lifetime of H-bonds in supercritical water is 0.2-0.5 ps [36,37]; see reviews by [38]).The total number of H-bonds increases as a function of density, consistent with the trends observed in Figures 4(a) and 5(a).Considering the water clusters linked by H-bonds, the simulations at a density of 0.02-0.025g/cm 3 yielded smaller clusters than the simulations conducted at higher density (0.03-0.09 g/cm 3 ).Overall, the  +  value is 1-1.5 at a density of 0.02-0.025g/cm 3 , and 2.5-3.5 at a density of 0.03-0.09g/cm 3 .
Figure 9 compares the hydration number derived from the MD simulations (RDF and geometrical H-bonds) with the average hydration number calculated from the (  * ))/ total , where   is the concentration of the hydrated CuCl species with  water molecules, n is the total number of hydrated CuCl species, and  total is the total concentration.thermodynamic properties for hydrated CuCl⋅nH 2 O 0 ( = 1-14) species derived by Migdisov et al. [3] from their solubility data.It can be seen that the hydration number derived from the MD simulation is significantly lower than that estimated from the data of Migdisov et al. [3].The reasons for this will be discussed in the next section.

Discussion
The overarching objective of this study was to compare the hydration behavior of the CuCl 0 complex derived from the MD with that determined experimentally by Migdisov et al. [3] for a similar temperature (350 ∘ C) and fluid density and to explain any differences in this behavior predicted by the two approaches.Hydration numbers for the [CuCl(H 2 O)] 0 moiety were derived from MD simulations by calculating the integral of the Cu-O paired distribution function and by counting the H-bonded water molecules using a geometrical definition of H-bonds (Table 2).As illustrated in Figure 9, the simulations confirm a trend of increasing hydration with increasing solution density, which tallies with the solubility measurements.However, the hydration numbers predicted using both methods are much lower in the simulations than those predicted by the experimentally based thermodynamic model.
Significantly, the MD simulations reveal that the H 2 O clusters are short-lived (a few picoseconds) and the second hydration shell is very dynamic.As a result, the hydration number varies greatly over short periods of time (Figures 3 and 7).This indicates that the hydration water in the second shell is only weakly bonded, even compared to similar simulations at higher density for CuCl and Au-Cl complexes [14,26].This raises an important question about the nature of the thermodynamic models that we use to model such solutions with density between those of (near) ideal gases and those of condensed solutions (density > ∼0.6 g/cm 3 ).
4.1.Nature of Hydrated CuCl 0 Species.Speciation-based thermodynamic models, as opposed, for example, to models that assume full dissociation of complexes (e.g., Pitzer, 1973), allow for successful extrapolations beyond experimental conditions because they include much important intrinsic information about the chemistry of the solutions [39].The concept of "speciation," however, can become ambiguous in some situations, in which case, it is necessary to reconsider the significance and validity of the assumptions linking physical reality and thermodynamic models.For example, Driesner et al. [33] illustrated the formation of a large variety of NaCl clusters in concentrated solutions.In such solutions, the speciation was defined by the time-averaged proportions of clusters such as NaCl 2 − and Na 2 Cl 2(aq) .Owing to the short residence times and variable geometry of different NaCl clusters, it is difficult to derive a speciation model with a precise geometry and stoichiometry for accurate extrapolation.
In the case of CuCl 0 in water vapor, our simulations reveal that the [H 2 O-CuCl] 0 moiety is very stable (residence time > 50 ps) with a well-defined linear structure, confirming the information obtained via solubility, in situ UV-Vis, and in situ XAS experiments at a density > 0.19 g/cm 3 [6-9, 12, 40-42].In contrast, the second hydration shell revealed in this study is highly disordered and dynamic (Figures 3 and 7).The hydration of the [H 2 O-CuCl] 0 moiety is, by nature, a statistical effect, and time-averaged hydration numbers only reflect the "averaged" nature of the interaction between the CuCl 0 cluster and the solvent molecules within certain time scales.

Hydrogen Bonding: Vapor versus High-Density Liquid.
To evaluate the effect of water-water interaction, we counted the total number of H-bonds in the solutions with different fluid density (or water fugacity).As shown in Figure 8(a), the running average of the number of H-bonds shows good convergence over the whole simulation time and increases with increasing fluid density (Table 2).The short lifetime of the [CuCl(H 2 O)  ] 0 clusters in low-density fluids reflects the weaker interaction between water molecules in the solvent at low-density compared to high-density, because changing configurations in low-density fluids requires breaking fewer H-bonds and hence a lower activation energy in low-density than in high-density fluids [43].
Following the approach of Kalinichev [43], the total Hbonds in the simulation box were expressed as "H-bond total per water molecule" by dividing total H-bonds with 55, the total number of water molecules in the simulation box (Table 2).Kalinichev [43] (using Monte Carlo and classical MD data presented in [38]) calculated the number of Hbonds per molecule in bulk water for a wide range in density over a temperature range of 25-1000 ∘ C. Their calculations showed a near-linear increase in the number of H-bonds as a function of density for bulk water, and a small decrease in their number with increasing temperature for a given density.The degree of H-bonding in the bulk fluid is an important parameter to consider when using a geometrical approach to define the hydration of a cluster; in high-density solutions, all water molecules interact via H-bonding, so that hydration defined in this manner becomes infinite.This happens when the number of H-bond per water molecule ∼1 at a density of ∼0.4 g/cm 3 over the 300-600 ∘ C temperature range [43], that is, slightly above the critical density of water (0.32 g/cm 3 ).At fluids density of ∼1 g/cm 3 , the theoretical percolation threshold for the network of hydrogen bonds in water is known to be ∼1.4-1.7, which has since been confirmed by multiple molecular simulations [44][45][46] and experiments [47,48].Our ab initio calculations at very low-density range (0.02-0.09 g/cm 3 at 350 ∘ C) indicated 0.09-0.25 H-bonds per water molecule.These values are consistent with the trends reported by Kalinichev [43].

Defining Microscopic Hydration Numbers and Comparison with Experiments.
Owing to the short lifetime and lack of definite geometry for second-shell complexes in low-density vapor, defining individual hydrated complexes [CuCl(H 2 O)⋅(H 2 O)  ] 0 is not a trivial task and will depend to some degree on the definition of what constitutes a cluster of particular stoichiometry (m).In this study, we used two approaches to quantify the hydration of the [H 2 O-CuCl] 0 cluster (Figure 1).(i) In the first approach (RDF approach), the hydration shell is obtained by calculating the number of water molecules in a defined sphere of fixed radius (5 Å) around the Cu(I) atom, based on the RDF.(ii) The second approach (Hbond approach) is to count the number of water molecules in the whole simulation box that are linked to the [H 2 O-CuCl] 0 cluster by H-bonds as defined by the criteria of Chialvo et al. [31] (dist D-H < 3.5 Å and angle D-H-A < 70 ∘ ).
Overall, these two approaches produce very similar trends and absolute values for the change of hydration as function of fluids density; however, both yielded smaller values than the numbers inferred from the experimental data (Figure 9).The numbers of water molecules contributing to the hydration of the [H 2 O-CuCl] 0 cluster calculated from the RDF and H-bond approach are similar, but the actual H 2 O molecules may differ, as illustrated in Figure 7.The reason for this is that the H-bond approach also considers the contribution of H 2 O molecules beyond 5 Å linked to the [H 2 O-CuCl] 0 cluster by H-bonds (Figures 7(a) and 7(d)).However, in some cases, the H-bond approach ignores the potential contribution of water molecules that are bonded to the Cu(I) ion via Van der Waals forces.For example, in Figure 7(b), some water molecules are linked to Cu by Van der Waals forces (grey dashed lines); these water molecules may be part of larger hydrated cluster and not considered by the H-bond approach.Consequently, both approaches provide a lower estimate of the actual hydration number (see four examples in Figure 7).Note that the combined number, n, of water molecules around [CuCl(H 2 O)  ] 0 in snapshots can total 8 (Figure 7), and this value is also very dynamic.
Although our analysis of MD data underestimates the hydration number of the CuCl complex predicted by the experimental solubility study [3], the MD results, nonetheless, show the same trend as those estimated from the experimental study.There are two likely reasons for the discrepancy between the hydration numbers derived from our ab initio MD simulations and the solubility experiments.One of these reasons, given the nature of ab initio MD (limited in simulation time and number of molecules), is that our calculations may have underestimated the formation of large, fairly stable clusters.The limitation in simulation size of ab initio MD means that our model does not reproduce the spatial heterogeneity in density, and thus it is difficult to evaluate all possible long-range effects arising from interactions of the CuCl 0 complex with solvent molecules having a large radius.For example, in the solubility study, the results of experiments at the highest density predicted a hydration number of 14 [3], whereas the MD analysis predicted a number of 4. This could be because a species such as [CuCl⋅14H 2 O] 0 does not possess a rigid stoichiometry and geometry and instead may represent an average hydration of state of CuCl 0 in a solution with molecules having a large radius (i.e., outside the hydration sphere of 5 Å calculated by RDF approach, Figure 7).This is undefined in MD analysis, as it represents a solution that is intermediate between a gas (mostly noninteracting molecules) and a solution (mostly interacting water molecules).Mixed classical and ab initio MD simulations, ground-proofed using the current work, may help to provide further insight into the nature of hydration and the longrange effect of the solvent in low-density fluids at these intermediate conditions.Another reason for the discrepancy between the hydration number predicted by our simulations and that predicted from solubility experiments is that the latter is sensitive to the choice of fugacity model used.In treating the experimental solubility data, Migdisov et al. [3] employed the Lewis-Randall model for gas mixtures (ideal mixture of real gases), in which the behavior of Cu-bearing clusters was assumed to obey the Ideal Gas Law.The authors assumed that most of the nonideality in the behavior of these species can be addressed by accounting for the formation of hydrogen bonds between gaseous CuCl and gas-solvent, H 2 O (cluster formation).The validity of this assumption, however, has not been verified, and we do not exclude the possibility that, in addition to hydration, other processes may contribute to the deviation of the formed moieties from ideal behavior.

Conclusions
The ab initio molecular dynamics (MD) simulations show that in a HCl-H 2 O vapor with density between 0.02 and 0.09 g/cm 3 at 350 ∘ C, Cu(I) is bonded with one water and one chloride ligand as a linear [H 2 O-Cu-Cl] 0 moiety, with a highly dynamic second-shell hydration shell that has short life-span configurations.The changes of hydration number as a function of water fugacity are qualitatively consistent with previous solubility under similar conditions, with lower hydration numbers than the experimentally inferred values [3].The discrepancy could be due to (i) our MD analysis underestimating the "thermodynamic" hydration number (e.g., ignoring the effects of density heterogeneity) in high temperature fluids; and (ii) uncertainty in the estimation of the fugacity coefficients for these highly nonideal "vapors" in the experiments.The calculated number of hydrogen bond per water molecule decreases with decreasing water density, indicating weak interaction between water molecules in such low-density fluids.This study provides the first theoretical confirmation that beyond-first-shell hydrated metal complexes play an important role in metal transport in lowdensity hydrothermal fluids, even if it is highly disordered and dynamic in nature.

Figure 1 :
Figure 1: A snapshot of the CuCl(H 2 O) 0 complex and the surrounding water.

Figure 2 :
Figure 2: The Cu-O radial distribution function (RDF, solid lines, left axes) and the coordination number (CN, dashed lines, right axes) versus the Cu-O distance.(a) The RDF of Cu with all O atoms in the simulation.(b) The RDF of Cu with O atoms in the solvent water (excluding the water molecule in the [CuCl(H 2 O)] clusters, second-shell hydration).

Figure 3 :
Figure 3: Dynamic hydration number of the CuCl 0 complex (within 5 Å of Cu) as a function of simulation time at different values of density (the density values are indicated at the top right in each diagram).
2 O)] 0 Complex.The radial distribution functions of Cu-O at different densities for all simulations times are plotted in Figure 2 (left -axes), together with their integrals representing the coordination number (right -axes).

Figure 2 (
a) shows the RDF of Cu with all O in the simulation box.The first peaks in Figure 2(a) indicate Cu-O bond distances of 1.91-1.93Å and correspond to the H 2 O bonded directly to Cu(I).There are also weak peaks at ∼3.5 Å that indicate the existence of a loosely bonded second water shell around the [CuCl(H 2 O)] 0 complex.This peak is highlighted in Figure 2(b), which shows the RDF of Cu with O calculated after excluding the first-shell water molecules.The hydration numbers of Cu at different densities are listed in Table Figure 4: H-O radial distribution function (RDF, solid lines, left axes) and the coordination number (CN, dashed lines, right axes) versus the H-O distance.(a) shows the RDF of all H with all O atoms in the simulation; (b) shows the RDF of H in the [CuCl(H 2 O)] 0 cluster with O in the solvent water (excluding the water molecule in the [CuCl(H 2 O)] 0 cluster).
Figure 5: O-O radial distribution function (RDF, solid lines, left axes) and the coordination number (CN, dash lines, right axes) versus the O-O distance.(a) All O in the simulation box; (b) O in the binding water (in [CuCl(H 2 O)] 0 clusters) with O in the solvent water (excluding the water molecule in the [CuCl(H 2 O)] 0 cluster).

Figure 6 :
Figure 6: Radial distribution function (RDF, solid lines, left axes) and the coordination number (CN, dash lines, right axes) as a function of the Cl-O distance of Cl versus O in solvent water (excluded water molecule in [CuCl(H 2 O)] 0 cluster).

Figure 7 :
Figure 7: Hydrogen bonds (red dashed lines) and the sphere of O within 5 Å of Cu (shown by the grey dashed lines) in the simulation for a density of 0.08 g/cm 3 .(H 2 O) within 5 Å of Cu is the hydration number calculated from the RDF and its integral; H-bonds were picked up when the D-H distance was less than 3.5 Å and the D-H-A angle less than 70 ∘ ; H-bonds total is the total number of H-bonds in the solution; (H 2 O) in an H-bond cluster is the number of water molecules in the [CuCl(H 2 O)  ] 0 cluster linked by H-bonds; the combined (H 2 O) value counts all the water molecules within 5 Å of Cu and water molecules beyond 5 Å that are linked by H-bonds.

GeofluidsFigure 8 :Figure 9 :
Figure 8: Running average for the total number of H-bonds (a) and number of water molecules surrounding the [CuCl(H 2 O)] 0 cluster linked by H-bonds (b) at different density.

Table 1 :
Simulation details for a model at 350 ∘ C with 1 Cu, 1 Cl, and 55 H 2 O in the box, and the geometry of the first-shell CuCl(H 2 O) complex.

Table 2 :
Number of H 2 O molecules and H-bonds from MD simulations for different density and water fugacity at 350 ∘ C.