Ab initio molecular dynamics free energy study of enhanced copper (II) dimerization on mineral surfaces

Understanding the adsorption of isolated metal cations from water on to mineral surfaces is critical for toxic waste retention and cleanup in the environment. Heterogeneous nucleation of metal oxyhydroxides and other minerals on material surfaces is key to crystal growth and dissolution. The link connecting these two areas, namely cation dimerization and polymerization, is far less understood. In this work we apply ab initio molecular dynamics calculations to examine the coordination structure of hydroxide-bridged Cu(II) dimers, and the free energy changes associated with Cu(II) dimerization on silica surfaces. The dimer dissociation pathway involves sequential breaking of two Cu2+-OH− bonds, yielding three local minima in the free energy profiles associated with 0-2 OH− bridges between the metal cations, and requires the design of a (to our knowledge) novel reaction coordinate for the simulations. Cu(II) adsorbed on silica surfaces are found to exhibit stronger tendency towards dimerization than when residing in water. Cluster-plus-implicit-solvent methods yield incorrect trends if OH− hydration is not correctly depicted. The predicted free energy landscapes are consistent with fast equilibrium times (seconds) among adsorbed structures, and favor Cu2+ dimer formation on silica surfaces over monomer adsorption.


S4. SUPPLEMENTARY NOTES 4
This section discusses the extraction of ∆G dimer from ∆W (R). ∆G dimer is related to ∆W (R) via the expression where the one-dimensional integral is over an appropriate 3-dimensional, R-dependent volume element Ω(R) dR, and V ref is the spatial volume associated with the standard state (e.g.,

S5
that of a solute at 1.0 M concentration, which is partly associated with the translational entropy of cation desorption or dimer dissociation).
Details of the AIMD trajectories used to calculate the potential of mean force are given in Table S1. The total trajectory lengths associated with the three configurations in Fig. 1a-c of the main text, excluding equilibration times, are 357 ps, 302 ps, and 456 ps, respectively.
In the Cu-dimer in pure water simulations (Fig. 1c of  For our dimerization problem, which has a highly non-trivial reaction coordinate R (Fig. 2 of the main text), Ω(R) is difficult to construct rigorously. As shown in Fig. S2, in the small R region of Fig. 3 of the main text, R is almost independent of R Cu−Cu . Our attempt to reconstruct a ∆W (R Cu−Cu ) from ∆W (R) fails because there is almost no volume element associated with R Cu−Cu in the initial stage of reaction. In contrast, in the asymptotic region where the two Cu 2+ are well-separated, R Cu−Cu is the correct coordinate to use; it integrates and t tot is the trajectory length used in sampling statistics in picosecond. δ, in meV, is the estimated uncertainty (one standard deviation) in each window.
In the case of Fig. 1a

S7
to a ∆W (R)-adjusted volume that is to be divided by V ref to yield a dimensionless quantity inside the logarithm expression (Eq. S1).
To accommodate these disparate requirements, we use a two-step procedure. (1) We perform one-dimensional integration of Eq. S1 over the inner (smallest R) and outer (largest R) free energy basins. Since the integrand in Eq. S1 exponentially decreases with ∆W (R) as R is varied from its local minimum, the predictions are insensitive to the precise 1-  bonds, but these k B T log(2) or k B T log(4) factors are small and are neglected. Recall that the main focus of this work is to demonstrate that Cu 2+ dimerization preferentially occurs on silica surfaces relative to liquid water. We expect some cancellation of uncertainties when comparing AIMD predictions for the three cases in the main text using the same Ω(R) definition, and the overall effect on the relative ∆G dimer with and without silica would benefit from cancellation of uncertainties. The overall systematic uncertainty in ∆G dimer in this approach is likely less than 0.1 eV, which is smaller than the discrepancy between the AIMD ∆G dimer reported in the main text and the corresponding g09 predictions.

S5. SUPPLEMENTARY NOTES 5
This section contains further discussions of the reaction coordinate R. We chose reaction coordinate R, further illustrated in Fig. S4, after rejecting others via trial-and-error, or by examining the energetics associated with adding a water to the Cu hydration shells. An analogy can be made with the early days of AIMD PMF modeling of acid-base reactions. The pKw of liquid water was first computed using a distance constraint. 6 A coordinate constraint approach 7 improved upon it. One of the present authors applied a 4-body coordinate for surfaces. Subsequently, Sprik and coworkers developed and perfected a statistical mechanical approach that annihilates the acid proton from the vicinity of the acid molecule and moves it far into the aqueous region. 8 Thus the later "reaction coordinates" improved upon the earlier ones; without the earlier work the development of the subsequent, improved coordinates might have been much delayed. It is in this spirit that we believe our present work is a useful, pioneering, if imperfect, contribution to computing the free energy of metal cation dimers. to an apparent kink in the overall PMF curve, we add a window with R o in between those two, with a larger A o , and rerun. Table S1 shows that the Fig. 3b calculation has a few of these added intermediate window. This is because we initially chose R o in two successive windows which are too far apart and the overlapping bin has small probabilities.
To illustrate some of the choices involved, the inset of Fig. S6 shows that two initial choices of adjacent windows exhibit discontinuous slopes in P (R). The two windows in question had of freedom that requires a secondary (quasi-2D) umbrella sampling.
In summary, we manually but carefully work match the P (R)s from different windows.

S6. SUPPLEMENTARY NOTES 6
This section discusses the reversibility of the reaction coordinate R. The ∆W (R) predictions depicted in Fig. 3 of the main text are obtained by starting from the compact dimer configuration (R Cu−Cu ∼2.9Å) and incrementally increasing the constraint R o value to dissociate the dimer. Here we show that this procedure is piecewise reversible. This demonstration is particularly important for our complex reaction coordinate R, which needs to specify tagged OH − /H 2 O; it is therefore vulnerable to those ions or molecules diffusing S12 away from the reaction zone, rendering R irrelevant. and R o =1.7Å windows, respectively. The results are similar. We note that, when we restart from a R o =2.0Å window configuration (two R o windows away) and attempt to recover the R o =1.4Å window behavior, that effort fails to recover the original R o =1.4Å window ∆W (R) because one of the bridging OH − groups that makes up the original R coordinate has acquired a H + from the surrounding water, and has diffused away from the bridging position. To recover the ∆W (R) statistics, we would need to designate a different H 2 O or OH − group as part of R. This suggests that our PMF calculation with coordinate R is piecewise reversible, but care must be exercised when assuming that it is globally reversible. Note that such concerns are not limited to the umbrella sampling method used herein; related approaches like metadynamics face similar reversibility issues.
In the case of the Cu-dimer horizontally adsorbed on silica surfaces (Fig. 1a/3a of the main text), another complication arises. As mentioned in the main text, the initial configuation has each Cu 2+ coordinated to a SiO − surface group. With sufficient equilibration in the smallest R o windows, one of the Cu 2+ becomes detached from the surface, and is coordinated to H 2 O and OH − only. As the Cu 2+ dimer starts to dissociate (R>∼0Å), however, it becomes favorable for both Cu 2+ to be coordinated to the surface. The AIMD/PMF trajectory lengths used are not always sufficient to "equilibrate" these two configurations (i.e., to allow the multiple occurances of reversible coordination of one of the Cu 2+ to the surface). Such complications are not infrequently encountered when dealing with the complex energy landscape in AIMD/PMF simulations of water-material interfaces. Fortunately, we find that ∆W (R) does not strongly depend on whether one or both Cu 2+ are bound to surface SiO − groups. Fig. S9a-b depict snapshots of two trajectories at R o =-0.50Å, initiated from R o windows which are smaller and larger, respectively. Throughout their respective AIMD trajectories, one and two Cu 2+ are coordinated to the silica surface. In other words, full equilibration of the Cu 2+ configuration has not been achieved. Fortunately, Fig. S8b shows that the ∆W (R) along these two trajectories are almost identical. Therefore ∆W (R) is not strongly affected whether one or both Cu 2+ is/are initially coordinated to the silica surface. to the two snapshots in Fig. S9a-b, respectively.

S7. SUPPLEMENTARY NOTES 7
This section provides more g09 calculation details. We consider general dimerization "reactions" of the form where all H 2 O and OH − species are understood to be in the Cu 2+ first hydration shells. In In the main text, we state that the implicit solvation of OH − species is the reason the cluster approach (Fig. 4  Explicit water-explicit water interactions are also well-represented using the cluster-plus-PCM approach. Thus we next consider where the left side represents two isolated H 2 O molecules hydrated by PCM and the right is a water dimer hydrated by PCM. ∆G=-0.03 eV after accounting for entropic corrections S15 discussed above. This near-zero value again shows that PCM is a good approximation of water-water hydrogen bonding.
In contrast, the bare OH − ion is not accurately "solvated" by PCM. Consider where all species are PCM-solvated and the cluster on the right side is depicted in Fig. S10.
After making the entropic corrections discussed above, the free energy change of the above "reaction" is +0.42 eV, which is significantly different from the ideal 0.0 eV value. This suggests that PCM over-stabilizes the bare OH − anion. The different configurations of the Cu 2+ dimer in Fig. 4  Instead, the coordination number of each OH − can be deduced using AIMD simulations.
As mentioned in the main text, the PCM cavity size can also be optimized to yield more accurate results. 5 This manuscript is focused on applying AIMD/PMF methods; we leave such development of cluster calculations to future work. S16

S8. SUPPLEMENTARY NOTES 8
This section compares the DFT method used in the main text with more accurate methods. Table S2 lists the Cartesian coordinates of "compound 1" 9 predicted using our triplet states DFT/PBE calculations, with the Gaussian suite of programs and a lanl2dz/6-311+G(d,p) basis set. This molecule is studied in Singh et al. 9 using more accurate, TPSS functional. Table S3 compares the bond lengths and bond angles predicted using these methods. DFT/PBE overestimates all bond lengths, which it is known to do. It also misses the symmetry-breaking in the bonding environment of the two Cu atoms. We believe the more pertinent issue is the energetics of Cu dimer separation in compound 1; such energy comparisons should be conducted in the future.

S9. SUPPLEMENTARY NOTES 9
This section documents the optimized cluster coordinates. S17