Isomerism of Trimeric Aluminum Complexes in Aqueous Environments

Article pubs.acs.org/JPCB

Isomerism of Trimeric Aluminum Complexes in Aqueous Environments: Exploration via DFT-Based Metadynamics Simulation Giorgio Lanzani,†,# Ari P. Seitsonen,*,‡,§ Marcella Iannuzzi,‡ Kari Laasonen,*,∥ and Simo O. Pehkonen*,⊥ †

Thule-Institute, University of Oulu, P.O. Box 7300, FI-90014 University of Oulu, Finland Institut für Chemie, University of Zurich, Winterthurerstrasse 190, CH-8057 Zurich, Switzerland § Département de Chimie, École Normale Supérieure, 24 rue Lhomond, F-75005 Paris, France ∥ Department of Chemistry, Aalto University, P.O. Box 16100, FI-00076 Aalto, Finland ⊥ Department of Environmental and Biosciences, University of Eastern Finland, PL 1627, FI-70211 Kuopio, Finland ‡

S Supporting Information *

ABSTRACT: The chemistry of aluminum or oxo-aluminum in water is still relatively unknown, although it is the basis for many chemical and industrial processes, including flocculation in water treatment plants. Trimeric species have a predominant role in the formation of the Keggin cations, which are the basic building blocks of aluminum-based chemicals. Despite this, details of the structural evolution of these small solvated clusters and how this is related to the processes leading to the formation of larger aggregates are still an open issue. To address these questions, here, we have applied the metadynamics (MTD) simulation technique [Barducci, A.; et al. Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2010, 1, 826−843] with density functional theory-based molecular dynamics to disclose the dynamics and structural conversions of trimeric aluminum complexes in an aqueous environment. The existence of a variety of competing metastable conformations, for example, book-like, cyclic boat, and linear shape conformations, is revealed in the MTD simulation. Furthermore, equilibrium simulations of the various intermediate states encountered along the MTD trajectory are used to assess their (meta)stability, determine the rearrangement of the OH ligands, and discuss the role of the solvating water.

■

INTRODUCTION Aluminum is one of the most widely used metals and one of the most common metals in the Earth’s crust. It plays a very important role in many chemical and industrial processes, including alumina-based sol−gel synthesis and hydrolysis, and it is used in large quantities as a flocculant chemical in water treatment plants all over the world. Because of its natural abundance, aluminum is commonly regarded as a harmless element, yet high levels of the water-soluble forms of aluminum may negatively affect terrestrial and aqueous life in different ways.1 When the level of Al(III) is kept under control, its hydrolyzed species can be widely and safely used for water purification purposes.2 Despite this, there are still several factors concerning the structures and the chemical composition during coagulation that need to be clarified to maximize the efficiency of these chemicals in an environmentally safe use. In particular, any improvements in the speciation of various aluminum structures and the understanding of the properties of the hydrolyzed Al species formed in water would be a significant contribution in understanding the aqueous oxo-aluminum chemistry. The relative stability of different polyaluminum complexes is highly sensitive to the Al(III) concentration and pH. The influence of these variables on the speciation of Al compounds © XXXX American Chemical Society

in an aqueous environment has been previously studied by means of electrospray ionization mass spectrometry (ESIMS),3−5 solid-state Al-nuclear magnetic resonance (NMR) spectroscopy,6,7 wet scanning electron microscopy (WSEM), tapping-mode atomic force microscopy,8 and computational chemistry techniques.9−12 Using ESI-MS techniques, Sarpola et al. have observed several different anionic and cationic aluminum complexes ranging from monomeric to Al30 species,3,4,11 yet the trimeric species were found to be among the most stable,13 which is in good agreement with the earlier ESI-MS results reported in the literature.14−17 It has been suggested that the core structure of these clusters would be a compact semispherical one, where three aluminum atoms are connected with OH bridges forming a six-membered ring with one common central oxygen O ligand.9 Four of such ring systems can be identified as a substructure of the δ isomer of the Keggin cations.18,19 Despite the multiple techniques used for studying the aquatic aluminum chemistry, many phenomena could not be unambiguously resolved yet. In fact, it is well known that the Received: August 11, 2016 Revised: October 19, 2016 Published: October 21, 2016 A

DOI: 10.1021/acs.jpcb.6b08112 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Figure 1. Geometry used for modeling the (a) cyclic and (b) linear isomers of the aqueous Al trimer (Al, pink; O, red; and H, white). Ball-and-stick model and lines are used to display the atoms in the cluster and in the solvating water, respectively. Greek letters denote the Al and O atoms involved in the cyclization reaction.

cluster cation (Al13O4(OH)24(H2O)12)7+, which acts as a key player in the flocculation process, presents a δ Keggin structure with a tetrahedral Al atom at its center, coordinated via four oxygen atoms to as many trimeric Al clusters. The formula can be thus expressed more fittingly as (AlO4Al12(OH)24(H2O)12)7+20 and it is generally called the Al13 ion. To understand the role of the trimer species in the formation of the δ Keggin structure, it is essential to confirm their stability in aqueous environments and to understand whether these species prefer to assume a linear (no Al−O−Al− O−Al−O ring present) or cyclic conformation. In this study, we have performed metadynamics (MTD)21 simulation on the basis of ab initio molecular dynamics (AIMD)22 of fully solvated trimeric aluminum complexes to observe their chemistry in water. This approach provides detailed information on the structures and structural transitions of the trimer without explicitly imposing a reaction path to follow, as was done in our previous work on the basis of the constrained dynamics approach.10,23 Particular attention has been devoted to the characterization of the linear and cyclic structures, as shown in Figure 1.

We have performed AIMD simulations on the basis of density functional theory (DFT) for the calculation of the electronic structure, providing the forces on the ionic cores needed to perform the MD. The starting structures of the trimeric Al complex have been first optimized in the gas phase and then inserted into a cavity in a cubic cell containing 141 water molecules and equilibrated at a finite temperature. Periodic boundary conditions with a lattice constant of 16 Å are applied throughout. The simulation cell is too small, or, contains too few water molecules to allow us to address macroscopic properties, such as pH or concentration, because the concentration can be addressed only in steps of 0.393 mol/ L. On the other hand, the cell is sufficiently large to isolate the cluster from its periodic images and investigate its dynamics in solution. The largest extent of the Al−Al distance is 8 Å and the smallest Al−Al distance over the periodic images is about 9 Å. Even taking into account the OH− groups, there are no significant artificial interactions between the periodic images and therefore it is highly unlikely that these would cause bias to the results shown here. The atomic system is in the charge state +1 and a homogeneous background charge neutralizes the simulation cell. All simulations are carried out using the Quickstep module26 of the CP2K program.27,28 It relies on the pseudo-potential approximation and Gaussian-plane wave (GPW) dual basis set methodology:29 Gaussian basis sets of triple-valence quality (TZVP-GTH)30 have been selected for the representation of the molecular orbitals. For the plane wave (PW) expansion of the charge density, a cutoff of 280 Ry is used. The generalized gradient-corrected Becke−Lee−Yang−Parr density functional is used for the calculation of the exchange and correlation terms.31,32 This functional has been already proven to be accurate in describing the structural characteristics and stability of aluminum(III) in aquatic environments.23 We note that this choice of the cutoff energy and exchange−correlation functional, without taking into account the London dispersion forces, has been shown to limit the accuracy of simulations in water.33 Here, we are, however, mainly interested in the local environment of the Al ions, which is not directly affected by these choices of simulation parameters. The structures of Al−O clusters have been earlier compared with different computational methods and were found to be remarkably similar.11

■

METHODS AND COMPUTATIONAL DETAILS Previous computational studies have already provided insight into the structures, stability, and electronic properties of aluminum compounds in an aqueous environment.10,23,24 The MTD approach employed in the present study is a methodology on the basis of molecular dynamics (MD) simulations and is used here to enhance the sampling of the phase space and eventually describe the mechanism of slow processes as a function of few collective variables (CVs).25 The selected CVs are assumed to provide a correct coarse-grained description of the processes of interest. In the space of the CVs, a historydependent potential is progressively constructed along the simulation to induce the system to visit still unexplored parts of the phase space by reducing the probability of coming back to already visited configurations. Such an MTD potential is obtained from the sum of small Gaussian contributions spawned along the MTD trajectory at regular time intervals. The size of the Gaussian beads depends on the characteristic fluctuations of the CVs and on the topography of the underlying free energy surface. B


Article


Figure 2. Values of the DD and RMSD CVs monitored along 24 and 22 ps of equilibrium AIMD-NVT simulations starting from the cyclic (left panels) or the linear cluster structure (right panels).

processes can become important for the stabilization of intermediate structures during the conformational changes of the trimeric complex. These are collective phenomena that cannot be efficiently described considering the coordination of single O atoms with H atoms. Therefore, we prefer to monitor the total number of zerofold (O), onefold (OH−), twofold (H2O), and threefold (H3O+) assemblies in the entire simulation cell. We have evaluated the number of hydrogen atoms associated with an oxygen atom by rounding off the O− H coordination number to the nearest integer. To guarantee an effective application of MTD, the CVs should differentiate relevant states along the transformations (i.e., equilibrium, intermediate, and transition states). A careful analysis of nonaccelerated MD trajectories, starting from the two known equilibrium structures of the solvated cluster, linear and cyclic, is useful to determine the CVs that properly distinguish different states in the phase space. We selected two CVs for the MTD simulation. One is the difference of distances (DD) between Al atoms defined as

All the MD simulations are carried out in a canonical, NVT ensemble at 300 K. To ensure fast temperature control, the system was coupled to a chain of Nosé−Hoover thermostats with a characteristic frequency of 2500 cm−1, which couples to the O−D stretching vibrations located around this frequency.34−36 The equations of motion are integrated using the velocity Verlet algorithm with a time step of 1 fs. The length of the MTD trajectory is 400 ps.37 We emphasize here that the “time” in MTD simulations, denoted as τ, is not directly related to the physical, or, real-time dynamics and it was not until recently that an attempt to extract the real-time dynamics from MTD simulations has appeared.38 To monitor the conformational changes of the solvated clusters along MD and MTD trajectories, we have employed simple geometrical quantities, such as distances and average coordination numbers of one atomic species with respect to a second one. The analytical definition of the generalized coordination number in our implementation is 1 Cab = Na

Na

∑ ia = 1

⎡ N 1 − 1 ⎢ b ⎢∑ Nb ⎢ i = 1 1 − ⎣ b

( ) ⎤⎥⎥ ( ) ⎥⎦ riaib

p

r0 riaib q r0

DD = d(Al γ −Alα) − d(Al γ −Alε)

(2)

It fluctuates around 0 Å when the structure is cyclic, whereas its absolute value can increase to about 2.5 Å for the linear form. Figure 2a,b shows the time evolution of the DD along two NVT simulations at 300 K starting from a cyclic and a linear conformer, respectively. To ensure that the system does not explore areas of the phase space that are not relevant for the description of the cyclization process, we place a restraining potential on the DD during the MTD simulation. This potential acts as Gaussian repulsive walls at DD = −4.5 and +4.5 Å and guarantees that fragments originating from the decomposition of the trimer do not get too far from each other. The second CV is defined in terms of the root mean square deviation (RMSD) with respect to a reference geometry with the distances of the present coordinates from those aligned in the reference squared. As references, we take the optimized cyclic and linear structures and restrict the calculation of the

(1)

where a and b label the two atomic species (or groups of atoms) containing Na and Nb elements, respectively, riaib is the actual distance between the two atoms indexed ia and ib, r0 is a reference distance that, together with the exponents p and q, determine the decay from 1 to 0 of the argument of the internal sum. For the coordination of Al atoms (a) with O atoms (b), the parameters are set to p = 250, q = 500, and r0 = 2.6 Å; for the O−H coordination, r0 = 1.25 Å. The radii were derived from the Al−O and O−H radial distribution functions in our system, which are shown in Figure S1 of the Supporting Information (SI). The corresponding quantity for the Al−Al distances is displayed in Figure S2. Most of the hydrolysis reactions occur as rearrangements of the hydrogen atoms or as associative hydration reactions. Such C


Article


Figure 3. Evolution of the three Al−Al distances during the simulation. The vertical lines separate time domains with different average distances and thus different structures of the Al−O ring. Characteristic geometries of the cluster during the simulation are reported as well (Al, pink; O, red; and H, white).

RMSD variable also takes into account the position of the O atoms belonging to the cluster. The rearrangements of the solvent, instead, are not directly addressed by this choice of CVs. Nevertheless, we expect that by activating structural changes in the trimeric complex with a smooth and slowly growing MTD potential, the correct behavior of the solvent is going to be reproduced. The MTD potential is obtained by adding Gaussian hills of a height of 2.7 meV. The width of the hills has been adjusted on the amplitude of the fluctuations of the respective CVs, and the values 1.0 for the DD and 0.15 for the RMSD are used in the MTD simulation.

RMSD to the atoms constituting the Al−O chain indicated with Greek letters in Figure 1; RMSDcy is the deviation from the cyclic geometry, whereas RMSDlin is the deviation from the linear one. The CV is defined as the normalized difference of the two, that is RMSD =

RMSDcy − RMSDlin RMSDcy + RMSDlin

(3)

As this function is normalized, it takes values close to −1 and +1 when the configuration is similar to the cyclic or linear configuration, respectively. Along the simulation, we have monitored the three Al−Al distances to follow the sequence of structural changes and to identify all relevant intermediate states. We remark that the Al− Al distances are less general variables than the combination of DD and RMSD because different combinations of the three distances might result in equivalent structures. Moreover, the

■

RESULTS Initialization of Cyclic and Linear Configurations. The first set of simulations have been carried out to sample the two known23 equilibrium states of [Al3O4]+(aq), namely, the linear and cyclic isomers shown in Figure 1. As mentioned above, D


Article


Figure 4. Evolution of the (a) DD and (b) RMSD along the MTD trajectory. (c) Correlation of the DD and RMSD during the different MTD time windows along the MTD trajectory. The legends show the MTD time in picoseconds.

The simulation starts at the cyclic isomer A (cf. Figure 3), which is stable during the first 10 ps of the MTD simulation. After this period, we observe a transition of this initial cyclic isomer into a quite rigid book-like conformation B. Both DD and RMSD, shown in Figure 4a,b, experience a small, step-wise deviation from their initial values when this transition occurs. At τ = 27 ps, the distance Alγ−Alε suddenly rises up to 6 Å to eventually fall back to the original value after the first 100 ps. On the other hand, during this period, the other two Al−Al distances do not change substantially, fluctuating around 3 Å. The transition at τ = 27 ps occurs as one Al−O link opens and the complex transforms into a more flexible conformation, herein named C, where the six-membered ring is replaced by a four-membered Al2O2 ring bonded to a −OH−Alγ(OH)3 group, as observed in the snapshot of the cluster structure visualized in the upper part of Figure 3. The corresponding change in the two CVs is more evident in this case: the DD seems to equilibrate around a value of −2.5 Å but the amplitude of the fluctuations is larger than that in the conformers A and B. The RMSD starts fluctuating strongly, suggesting that the −OH−Al(OH)3 group can bend strongly from the rigid Al2O2 unit and that the encountered structures are rather divergent from both the cyclic and the linear reference structures. The steady increase in the value of RMSD toward 1 shows that the geometries visited at τ ≈ 100 ps are the closest to the reference linear conformer. The next geometry D for τ = 104−140 ps is again close to a cyclic structure, with all three Al−Al distances around 3 Å. The cyclic geometry is evidenced by DD remaining close to zero. Two of the Al atoms constitute an Al2O2 unit with both

these simulations are useful to determine which parameters distinguish the different states along the simulation. The results from the analysis of the MD trajectories (cf. Figure 2) confirm the thermodynamic stability of both conformers: during the 22−24 ps long simulations, the systems equilibrate quickly (on a time scale of a few picoseconds) and no variation of the Al−O−Al connectivity is observed. The cationic aluminum ions keep their fourfold coordination in the aqueous environment. This is in good agreement with the general understanding that in water Al atoms form polymeric cations, in which the unit formula is a Al3+ central atom tetrahedrally surrounded by four oxygen atoms ([AlO4]5−).39 The DD and RMSD CVs monitored along the equilibrium MD simulations are shown in Figure 2. The CVs clearly differentiate between the two isomers. Structural Transformations from MTD Simulation. The MTD simulation has been started from the equilibrated cyclic structure using as CVs the DD and the RMSD, as described above. The time evolution of the three Al−Al distances is shown in Figure 3. The plot is partitioned into MTD time windows, which are characterized by different combinations of values of the three distances, that is, different structural features of the trimeric complex. In particular, when all the three distances have the same average, the cluster is cyclic, whereas when there are two smaller and one larger value, the cluster is in the linear form; one Al−Al oscillating around 3 Å can be associated with a simple Al−O−Al bond topology, and a dissociation of the trimer has occurred when two distances simultaneously exceed a value of about 4 Å. E


Article


Figure 5. Coordination of Al by the OHn groups along the MTD simulation. The sum over the oxygen species is somewhat lifted vertically for better visibility. Black, O; red, OH−; green, H2O; and orange, sum.

bridging ligands being OH− groups. Occasionally, one of those O ions is connected to two, more often to all the three Al atoms. The following geometry E is a cage structure formed by the two Al atoms in the Al2O2 unit connected via a further O− Alγ−O chain. Two of the Al−Al distances are slightly prolonged due to a further ligand attaching to Alγ when structure E is formed. The DD falls to a value −0.7, and the amplitude of the fluctuations remains small. For a short time, the structure passes again to a cyclic structure F that is similar to geometry D, being characterized by similar Al−Al distances and value of DD. At τ = 174 ps, the Alε−Alα distance increases to larger values into structure G, corresponding to the opening of one of the two rings of structure E. This leads to the formation of intermediate structures that are quite unstable and short-lived. There is a hydrogen bond from the H of one H2O of Alε onto an OH on Alα. With the Al2O2 ring and the third Al bound to them via the oxygen of a single OH group, structure G is similar to structure C. These structures are, however, different from the CV point of view because the latter do not include permutations of the atoms. The DD variable assumes positive values in this phase of the simulation, whereas the RMSD still fluctuates over a wide range of values but shows a tendency toward linearization of the structure by moving to values closer

to 1. After geometry G, the system visits structure H, which is built of quite a linear Al−O−Al−O−Al chain. The linearity is apparent in the large value in one of the three Al−Al−Al angles (cf. Figure S3). There is a hydrogen bond from one H2O ligand of Alα to the OH on Alγ at the center of the chain. At τ = 222−254 ps, denoted as structure I, two of the Al−Al distances are mainly above 4 Å, indicating that the trimer attempts to dissociate. This is prevented by the repulsive wall on the DD CV. As discussed above, the dissociated configurations are not of interest to us in this study, but the corresponding relatively deep minimum in free energy, due to large configurational entropy contribution, has to be filled before the MTD trajectory can move further and new structures are again visited. This happens at τ = 254 ps, when structure J, similar to structures C and G, is reached. This converts soon further, as all the three distances return to values of 3−3.5 Å at τ = 266 ps. In the MTD time window between 266 and 312 ps, the MTD trajectory visits other cyclic conformers. Again, the DD fluctuates between 0 and 1 and the RMSD becomes closer to −1. This indicates a structural similarity with the initial cyclic complex. The characteristic structure identified along this part of the trajectory is the cyclic conformer structure K. It is similar to structure A but differs from the other previous cyclic structures D, E, and F in that one F


Article


Figure 6. Number of OHn species along the MTD simulation.

sufficient to span the whole range of RMSD. Other regions instead are more interconnected and there are more transitions between them. The perfectly linear cluster, which would be identified with a value of RMSD of 1, is never visited, the closest approaches being in structures C and K. The MTD time windows at 104−140, 140−171, and 266−312 ps are characterized by narrow variations in the value of DD but large ranges of the value of RMSD. The latter two of these are located symmetrically with respect to DD, and in the second time period, the value of DD fluctuates around 0. At τ = 27− 104, 174−212, and 222−266 ps, both CVs experience relatively large variations in a quite short MTD time region. Role of the Solvent. To better understand the structural changes taking place during the MTD simulation, the Al−OHn connectivity and the number of OHn groups have been monitored in terms of the coordination numbers CAl−O and the total number of O atoms coordinated with zero to three H atoms, respectively. The resulting quantities are plotted in Figures 5 and 6, correspondingly. Because of the MTD bias potential, the system explores a large variety of Al−O−Al bonding configurations involving substantial rearrangements and reordering of the metal atoms in the chain. The time evolution of CAl−O in Figure 5, computed for individual Al atoms, shows that each cationic aluminum is

link in structure A between the Al atoms is a pure O atom, whereas in the three other structures, all the linking oxygencontaining ligands were OH groups. The Alγ−Alα distance is consequently somewhat shorter than the other two Al−Al distances from τ ≈ 280 ps. At τ = 312 ps, the cyclic structure opens again, first into a linear structure L and soon after to a new geometry, indicated by the simultaneous increase of Alε−Alα and Alγ−Alα distances, and it evolves into the dissociated dimer plus monomer fragments exemplified by geometry M. The fragments remain at relatively short distance due to the presence of the potential wall that becomes active again at DD values larger than 4.5 Å. The effect of the wall potential is clearly visible in the plot of the DD variable in Figure 4a. The fragments remain also connected via hydrogen bonds between the OH and H2O ligands on them and the central hydrogen atom often has an ambiguous character between the two OH groups. The map of the region explored in the two-dimensional CV space during the simulation is plotted in Figure 4c. Different colors have been chosen to indicate different time frames along the MTD simulation. This plot demonstrates that both CVs span the full range of relevant values. The area corresponding to the cyclic conformer is clearly separated from the rest and is visited only briefly at the beginning. The first 100 ps are G


Article


C. The oxygen atom, which is bonded to all the three Al ions in B, loses its bond to Alγ and becomes bridging between Alε and Alγ and captures one hydrogen atom from the solvent. One water molecule from the solvent binds as ligand to Alγ with a subsequent dehydrogenation of this water ligand, whereby one of its hydrogens moves to the solvent. During the reorganization process, the overall number of oxygen atoms in the cluster increases from 10 to 11 and the number of hydrogen atoms from 12 to 14. A reverse process takes place at τ = 270 ps when the system evolves from the linear structure J to the cyclic one K by the association of one water molecule from the proton of the bridging OH− group in structure J via a water bridge to Alε and one of its OH− ligands. On the other hand, the formation of stable H3O+ species has not been observed as demonstrated by the plot in the panel “H3O+” of Figure 6. Only fleetingly an O atom is coordinated by three H atoms, as these complexes are short-lived, being stable for tens of femtoseconds, and only rarely close to a picosecond. MD at Intermediate Structures. To test whether the isomers observed during the MTD simulation are real (meta)stable structures, AIMD-NVT simulations are performed using a characteristic snapshot configuration in these geometries as the starting point. We excluded the structures where the trimer dissociates, and structure A as it was already considered as an initial structure and discussed above and structure F as it occurs only very briefly and is similar to structure D. Well-equilibrated results are obtained during a simulation time of ≈40 ps with the duration of the simulation and the average Al−Al distances shown in Table 1; the time

always at least fourfold coordinated in the aqueous environment. This behavior is in good agreement with the general trend of Al to form a tetrahedron with four oxygen atoms [AlO4]5−.39 Nevertheless, because of its small size and a relatively large charge, the cation can also become easily fivefold or even sixfold coordinated. Such events are indeed frequently observed along the MTD trajectory when the trimeric complex undergoes structural rearrangements. Sixfold octahedral coordination occurs in the trimeric phase with two of the ligands on the sixfold Al ion bridge to the other Al ions, for example, in structures G and K as well as when one of the Al atoms is dissociating, for example, in structure M. A slow but systematic trend in the coordination of the Al ions is the change from the mainly four OH ligands of Alα at the beginning of the simulation to four H2O, or water, ligands at about τ = 200 ps and finally to six water ligands toward τ = 380 ps when Alα dissociates from the dimeric complex of Alε and Alγ. On Alε and Alγ, the type of ligands varies along the MTD simulation but does not undergo such a general change. Ionic O2−, without hydrogens attached, appears as a ligand only fleetingly and practically always bridging two Al ions, and we do not observe H3O+ ligands binding to the Al ions. Most of the lowcoordinated aluminum species, previously identified as stable species in the gas phase (see, e.g., ref 40), are never observed along our trajectory. Thus, they are unstable in the aqueous environment. Some of the changes in the Al−OHn coordination occur upon isomerization reactions of the trimeric complex. The first of such collective rearrangement is observed when structure A is converted into the book-like structure B and the Al−OHn coordination numbers of Alα and Alε increase from 4 to 5. The former has two own OH ligands and shares two other OH groups with two other aluminum atoms. Moreover, one single oxygen is coordinated to all the three Al atoms. Another large reshuffling of the coordinations occurs when the trimer switches from structure F to G or from I to J, and at about τ = 270 ps. More information about the oxygen species coordinating to the Al atoms is shown in Figures S4−S6. From the snapshots and a visual inspection of the MTD trajectory, we have observed that there is a series of structures characterized by the Al2O2 rings, showing a relatively high stability as they are often encountered along the trajectory. This finding is in good agreement with previous results of Saukkoriipi et al.24 We further observe that it is the OH group that is (most of the time) the bridging element between the two Al atoms, with a bridging oxygen atom only in structures A, K, and L of the trimer and remaining in structure M in the dimer unit when the trimer dissociates. Also, the preference of the OH bridges is in agreement with previous theoretical studies.23 The number of oxygen and oxygen−hydrogen species along the MTD simulation is monitored and shown in Figure 6. As expected, most of the O atoms form water molecules. There are, however, large fluctuations, and from our analysis we count that 103 of the initial 141 water molecules remain intact throughout the MTD simulation. Some water dissociation processes are observed, implying variations in the number of OH− and H3O+ groups as well. In particular, one large rearrangement process occurs at 26 ps when the number of OH− species increases from 6 to 8 because of the shift of one H from a water ligand of Alε to the bare oxygen atom. This reaction is related to the opening of the book-like structure B of the cluster, which is transformed into the quasi-linear structure

Table 1. Equilibrium MD Simulations from Snapshots in the MTD Trajectory: Simulation Time (in ps) and the Average Al−Al Distances (in Å) from Equilibrium MD Simulations Started at a Snapshot of Each of the MTD Time Windows initial structure

tMD [ps]

B C D E G H K L M

42.6 42.6 39.8 38.8 38.7 45.6 38.7 42.1 38.8

d(Al−Al) [Å] 3.49 5.00 3.24 3.26 3.45 3.44 3.17 3.39 3.20

2.91 2.97 2.96 2.99 4.73 6.19 3.50 4.77 6.40

description 2.81 3.39 3.13 3.35 2.98 3.34 3.50 3.57 5.45

stable stable stable converts into D stable stable stable stable stable

evolution in the Al−Al distances is shown in Figures S7−S9. During these simulations, the initial geometries equilibrate in a few picoseconds when the solvent slightly rearranges around the cluster, but the Al−O−Al or Al−OH connectivities do not change in most of the simulations and thus there are no changes in the solvent molecules. The only exception is the simulation starting from structure E, where the isomer converts into structure D after about 16 ps by Alγ losing one of its ligands and one of the O ions of the Al2OH2 bridge binds again fleetingly to all the three Al ions. Except for the small change in the geometry to accommodate this structural conversion, no large distortions occur neither during the transition nor in the following 20 ps of the simulation. Thus, we find that the structures, excluding E, indeed correspond to (meta)stable structures, at least on the time scale of our simulations of tens of picoseconds. H


Article

The Journal of Physical Chemistry B Inapplicability of Implicit Solvent. During the analysis, it turned out that the free energies are difficult to estimate from the MTD simulation, but in principle one could try to use static quantum chemistry calculations with an implicit solvation model to estimate the free energy of characteristic geometries along the MTD trajectory. To test this idea, we performed a geometry optimization procedure from the geometries in Figure 3. We took an initial geometry from each of the structures A-M and did geometry optimizations using an implicit solvation model (see SI for the details of the calculations). It turned out that the geometries from MTD configurations were very unstable in these calculations: in all optimizations, the structure of the cluster changed substantially, and in several cases we obtained unphysical results, for example, a hydrogen atom moved from oxygen to aluminum or we observed a formation of H2. Even if the idea that the free energies can be computed from the snapshot geometries is appealing, the geometries are so far from those in the simulation with the explicit water solvent that such an approach does not work, as the optimized structures cannot be used to evaluate free energies. This conclusion is in line with our earlier study,24 where we have shown that the gas-phase- and solvation-model-optimized structures of Al5OnHm clusters are not stable in explicit molecular water. Our test calculations here clearly show that the reverse case is also true.

stress that the choice of the CVs is nontrivial, as noticed along the course of the present study. In principle, with appropriate CVs, the MTD methodology applied here can be used for several different clusters. The selection and mapping of the high-dimensional coordinate space into two to three CVs are, however, a difficult task, complicated by the possibility of dissociation and rearrangement of water molecules. Overall, we can say that the Al trimer does not have a single structure in water. Yet, some of the structures might well be only transient states and their stability and existence probability depend on the relative free energies of the isomers so that few or even only one structure may be dominating. With our specific choice of the CVs, we are, however, not able to make any statements about the free energies. This needs to be tackled very carefully as the mapping of the multidimensional free energy surface to few CVs is very challenging, both due to the local coordinations of the Al ions and the hydration reactions with the solvent, rendering the pH as an important input variable determining the thermodynamic stability of different isomers.41

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.6b08112. Information about the Al−Al−Al angular distribution; coordination of the oxygen ions on the Al ions; Al−O, O−H, and Al−H radial distribution functions; histogram of Al−Al distances; Al−Al distances in the equilibrium MD trajectories, testing the (meta)stability of the structures observed; details of the calculations with the implicit solvent model (PDF)

■

CONCLUSIONS From the MTD simulations addressing the isomerism of aluminum trimer cluster in aqueous solution, several interesting observations have been made: the MTD technique offers advantages over previously applied constraint methods23 in enhancing the sampling of the configuration space. This is a highly nontrivial task because the coordination of the aluminum ions and their ligands changes during the simulation, arising from Al−OHn bond-formation and -breaking, and the solvent molecules participate in these reactions. We note that two of the structures obtained have been previously observed in implicit-solvent COSMO calculations.9 Th e aut ho rs f o un d in t h eir st ud y t ha t [ Al 3 O(OH)0−3Cl1−3(H2O)0−3]+ cations prefer core structures with a very compact three-dimensional cage-like C3v symmetry that characterizes also structure A in this study. Moreover, the peculiar stability of structures characterized by the presence of Al2O2 rings suggested by Saukkoriipi et al.9 is supported by the MTD simulations (structures C, D, E, G, and J). Many of the novel (meta)stable structures, especially the cyclic clusters (D/F and K), have low symmetry structures resulting from reactions of water molecules with the Al complex. We did not observe the fast structure breaking observed in the case of aluminum pentamer found by Saukkoriipi and Laasonen,24 where the highly symmetric and compact gas-phase structures open in the aqueous phase. It can thus be that the best structural motifs of aluminum-oxo clusters are size dependent. On the other hand, also in experiments, the trimer Al3(aq) is found to be very stable and the trimeric unit is smaller than the pentamer; therefore, the latter is possibly only a transient intermediate toward the larger compact structures, such as the Keggin Al13 complex, whereas the trimer remains intact longer due to its stability. Also, many of the trimer structures are quite open, which is in agreement with the results of Saukkoriipi and Laasonen.24 Long MTD simulations were needed to map the CV configuration space of small oxo-aluminum clusters here. We

■

AUTHOR INFORMATION

Corresponding Authors

*E-mail: Ari.P.Seitsonen@iki.fi. Phone: +33 1 44 32 32 78 (A.P.S.). *E-mail: Kari.Laasonen@aalto.fi. Phone: +358 40 557 00 44 (K.L.). *E-mail: Simo.Pehkonen@uef.fi (S.O.P.). Present Address #

Sapa Precision Tubing A/S, Hydrovej 6, DK-6270 Tønder, Denmark (G.L.) Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We are thankful to J. Rämö, C.A. Morrison, and A. Sarpola for useful discussions. We thank the anonymous reviewer of our manuscript for suggesting the optimization calculations using implicit solvent. The work was carried out under the Academy of Finland, Project AquAlSi. We acknowledge the generous computational capacity provided by Centre for Scientific Computing (CSC, Finland) and Centro Svizzero di Calcolo Scientifico (CSCS, Ticino, Switzerland), which were awarded for the calculations. Traveling expenses were financed by the Research Council of the University of Oulu.

■

REFERENCES

(1) Schwoyer, W. L. K. Polyelectrolytes for Water and Wastewater Treatment; CRC Press: Boca Raton, FL, 2000. I


Article

The Journal of Physical Chemistry B (2) Tchobanoglous, G.; Burton, F. Wastewater Engineering; McGrawHill: Singapore, 1991. (3) Sarpola, A.; Hietapelto, V.; Jalonen, J.; Jokela, J.; Laitinen, R. S. Identification of the hydrolysis products of AlCl3·6H2O by electrospray ionization mass spectrometry. J. Mass Spectrom. 2004, 39, 423− 430. (4) Sarpola, A. T.; Hietapelto, V. K.; Jalonen, J. E.; Jokela, J.; Rämö, J. H. Comparison of hydrolysis products of AlCl3·6H2O in different concentrations by electrospray ionization time of flight mass spectrometer (ESI TOF MS). Int. J. Environ. Anal. Chem. 2006, 86, 1007−1018. (5) Urabe, T.; Tanaka, M. Time-resolved analysis of hydrolytic aluminum species in the formation of the tridecamer using electrospray ionization mass spectrometry. Rapid Commun. Mass Spectrom. 2011, 25, 2933−2942. (6) Huang, L.; Tang, H.-X.; Wang, D.-S.; Wang, S.-F.; Deng, Z.-J. Al(III) speciation distribution and transformation in high concentration PACl solutions. J. Environ. Sci. 2006, 18, 872−879. (7) Exley, C. Reflections upon and recent insight into the mechanism of formation of hydroxyaluminosilicates and the therapeutic potential of silicic acid. Coord. Chem. Rev. 2012, 256, 82−88. (8) Lin, J.-L.; Chin, C.-J. M.; Huang, C.; Pan, J. R.; Wang, D. Coagulation behavior of Al13 aggregates. Water Res. 2008, 42, 4281− 4290. (9) Saukkoriipi, J.; Sillanpäa,̈ A.; Laasonen, K. Computational studies of the cationic aluminium(chloro) hydroxides by quantum chemical ab initio methods. Phys. Chem. Chem. Phys. 2005, 7, 3785. (10) Saukkoriipi, J. J.; Laasonen, K. Density functional studies of the hydrolysis of aluminum (chloro)hydroxide in water with CPMD and COSMO. J. Phys. Chem. A 2008, 112, 10873−10880. (11) Sarpola, A. T.; Saukkoriipi, J. J.; Hietapelto, V. K.; Jalonen, J. E.; Jokela, J. T.; Joensuu, P. H.; Laasonen, K. E.; Rämö, J. H. Identification of hydrolysis products of AlCl3·6H2O in the presence of sulfate by electrospray ionization time-of-flight mass spectrometry and computational methods. Phys. Chem. Chem. Phys. 2007, 9, 377−388. (12) Coskuner, O. Single ion and dimerization studies of the Al(III) ion in aqueous solution. J. Phys. Chem. A 2010, 114, 10981−10987. (13) Rämö, J. H.; Sarpola, A. T.; Hellman, A. H.; Leiviskä, T. A.; Hietapelto, V. K.; Jokela, J. T.; Laitinen, R. S. Colloidal surfaces and oligomeric species generated by water treatment chemicals. Chem. Speciation Bioavailability 2008, 20, 13−22. (14) Baes, C.; Mesmer, R. The Hydrolysis of Cations; John Wiley and Sons: New York, 1976. (15) Ö hman, L.-O.; Forsling, W. Equilibrium and structural studies of silicon and aluminium in aqueous solution. 3. Potentiometric study of aluminium(III) hydrolysis and aluminium(III) hydroxo carbonates in 0.6 M Na(Cl). Acta Chem. Scand., Ser. A 1981, 795−802. (16) Ö hman, L.-O. Equilibrium and structural studies of silicon(IV) and aluminum(III) in aqueous solution. Inorg. Chem. 1988, 27, 2565− 2570. (17) Shock, E. L.; Sassani, D. C.; Willins, M.; Sverjenski, D. A. Inorganic species in geologic fluids: Correlations among standard molal thermodynamic properties of aqueous ions and hydroxide complexes. Geochim. Cosmochim. Acta 1997, 61, 907−950. (18) Casey, W. H. Large aqueous aluminum hydroxide molecules. Chem. Rev. 2006, 106, 1−16. (19) Hench, L.; West, J. Chemical Processing of Advanced Materials; John Wiley and Sons: New York, 1992. (20) Son, J. H.; Kwon, Y.-U.; Han, O. H. New ionic crystals of oppositely charged cluster ions and their characterization. Inorg. Chem. 2003, 42, 4153−4159. (21) Laio, A.; Parrinello, M. Escaping free-energy minima. Proc. Natl. Acad. Sci. U.S.A. 2002, 99, 12562−12566. (22) Hutter, J.; Marx, D. Ab Initio Molecular Dynamics; Cambridge University Press: Cambridge, U.K., 2009. (23) Lanzani, G.; Sarpola, A.; Saukkoriipi, J.; Laasonen, K.; Morrison, C. A.; Slater, B.; Rämö, J.; Pehkonen, S. O. Study of the stability of aluminium trimeric clusters in aqueous solutions. Mol. Simul. 2012, 38, 934−943.

(24) Saukkoriipi, J.; Laasonen, K. Theoretical study of the hydrolysis of pentameric aluminum complexes. J. Chem. Theory Comput. 2010, 6, 993−1007. (25) Barducci, A.; Bonomi, M.; Parrinello, M. Metadynamics. Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2011, 1, 826−843. (26) VandeVondele, J.; Krack, M.; Mohamed, F.; Parrinello, M.; Chassaing, T.; Hutter, J. Fast and accurate density functional calculations using a mixed Gaussian and plane waves approach. Comput. Phys. Commun. 2005, 167, 103−128. (27) Hutter, J.; Iannuzzi, M.; Schiffmann, F.; VandeVondele, J. cp2k:atomistic simulations of condensed matter systems. Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2013, 4, 15−25. (28) CP2k Developers Group under the Terms of the GNU General Public Licence (GPL), 2014. See http://www.cp2k.org/. (29) Lippert, G.; Hutter, J.; Parrinello, M. A hybrid Gaussian and plane wave density functional scheme. Mol. Phys. 1997, 92, 477−488. (30) VandeVondele, J.; Hutter, J. Gaussian basis sets for accurate calculations on molecular systems in gas and condensed phases. J. Chem. Phys. 2007, 127, No. 114105. (31) Becke, A. D. Density-functional exchange-energy approximation with correct asymptotic behavior. Phys. Rev. A 1988, 38, 3098. (32) Lee, C.; Yang, W.; Parr, R. Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density. Phys. Rev. B 1988, 37, 785. (33) Jonchiere, R.; Seitsonen, A. P.; Ferlat, G.; Saitta, A. M.; Vuilleumier, R. van der Waals effects in ab initio water at ambient and supercritical conditions. J. Chem. Phys. 2011, 135, No. 154503. (34) Nosé, S. A unified formulation of the constant temperature molecular dynamics methods. J. Chem. Phys. 1984, 81, 511. (35) Hoover, W. G. Canonical dynamics: Equilibrium phase-space distributions. Phys. Rev. A 1985, 31, 1695. (36) Martyna, G. Remarks on “Constant-temperature molecular dynamics with momentum conservation”. Phys. Rev. E 1994, 50, 3234. (37) We note that we continued out the simulation up to 800 ps but because no significantly new information was gained in the second half of the simulation, we only present here the first half. (38) Tiwary, P.; Parrinello, M. From metadynamics to dynamics. Phys. Rev. Lett. 2013, 111, No. 230602. (39) Shriver, D.; Atkins, P. Inorganic Chemistry, 3rd ed.; Oxford University Press: Oxford, U.K., 1999. (40) Demyk, K.; Van Heijnsbergen, D.; von Helden, G.; Meijer, G. Experimental study of gas phase titanium and aluminum oxide clusters. Astron. Astrophys. 2004, 420, 547−552. (41) Sarpola, A. The Hydrolysis of Aluminium, a Mass Spectrometric Study. Ph.D. Thesis, University of Oulu, 2007.

J


Isomerism of Trimeric Aluminum Complexes in Aqueous Environments

Recommend Documents