Solvation Structure and Dynamics of Aqueous Solutions of Au+ Ions: A Molecular Dynamics Simulation Study

Although gold (Au) as an element is precious and noble, its elemental as well as ionic form is of huge scientific importance in view of its applications in electrochemistry, nano-electronics and other related fields. We have studied structure and dynamics of aqueous solutions of gold ions (Au+) using molecular dynamics simulations. Using a modified LJ parameter set for the Au+ ions in water in our molecular dynamics simulations, we have established the hydration structure and dynamics of Au+ ions in terms of radial distributions, orientations and residence time of the nearest neighbours. Our results on peak position, height and coordination numbers are in much better agreement with those from the recent CPMD simulations. Relative orientation of the neighbours as obtained from the angular distributions suggests octahedral or trigonal bi-pyramidal structure of the solvation shell. Orientational distributions of dipoles and other molecular orientational vectors indicate that the hydrogen atoms of the water molecules are away from the central Au+ ion. The residence time calculated from the corresponding time correlation function is found to be reasonably high, indicating less exchange of water molecules between the first and second solvation shells of the Au+ ion. In essence, the present results are in much better agreement with the CPMD results as compared to other QM/MM and classical force-field-based simulations.


Introduction
From early civilisation, gold is one of the most precious metals due to its vast use in our life. The uses of gold in various areas include medicine [1,2], electronic appliances [3], optical devices [4][5][6], energy [7,8], water treatment [9][10][11] and so on. With increasing utilisation of gold in ever expanding electronics industries, accumulation of electronic wastes poses one of the biggest challenges to the environment. Therefore, recycling of the electronic waste and thereby recovering gold and other precious metals from it have been of utmost importance. In order to recover gold and other precious metals from the electronic wastes, dissolution of this metal in ionic form in an aqueous medium is necessary and therefore understanding the interaction of gold (and other precious metal) ions with water is extremely important. Previously, many researchers have investigated [12][13][14][15][16][17][18] gold in its neutral form namely in the form of nanoparticles, surfaces, clusters, etc. Classical force-field (FF)-based modelling [19][20][21] and simulation of gold in its ionic (Au + ) form are very rare. To the best of our knowledge, the first molecular dynamics (MD) simulation investigation of Au + ions in bulk water was reported by Rode and co-workers [22]. They have used hybrid quantum mechanics/molecular mechanics (QM/MM) simulations to study hydration behaviour of one Au + ion in its bulk aqueous solution. In this work, the first solvation shell of gold ion (extending up to 3.8 Å as obtained from their classical FF-based simulation) has been treated quantum mechanically, whereas usual force-field-based MM methodology has been used for the rest. The MM simulation employed in this study involves complicated potential terms with two-as well as three-body interactions. According to these hybrid simulation results, radial distribution function (RDF) of oxygen (OW) or hydrogen (HW) atoms of the water molecules around the Au + ion shows two solvation shells with the maxima at around 2.45 Å and 4.9 Å from the Au + ion. The distribution of nearest (1st shell) neighbour numbers predicts the Au + ion to be 3-, 4-, 5-, 6-and sevenfold coordinated with average coordination number of 4.7. The most probable coordination numbers obtained from this study are 4 and 5. The second shell of the Au + ion as obtained from this QM/MM simulation is quite large with an average population of 28 water molecules and the spread in the number of 2nd shell neighbours ranges from 26 to 37. The large spread in the coordination number of both the 1st (3 to 7) and 2nd (26 to 37) shells and the non-zero RDF values between the two peaks imply a very labile hydration structure for both the shells. As far as geometrical structure of the 1st shell is concerned, the QM/MM simulation has predicted a square planar structure. The octahedral geometry can be ruled out as the distribution of average O w -Au-O w angles which show a peak at 90º but not at 180º. As already mentioned, they have also identified rapid exchange of water molecules between two neighbouring shells. In order to improve their own results on hydration of Au + ion, the same group has reported [23] a study of aqueous solution of Au + ions using an improved quantum mechanical charge field molecular dynamics (QMCF) method. Owing to the observation of rapid exchange of water molecules between the 1st and the 2nd hydration shells in their previous study [22], the quantum region in this study was extended up to the 2nd solvation shell. However, the results of the later study are quite surprising. The Au + -Ow RDF shows an extra peak between the conventional 1st and 2nd peaks and, considering the nearest 1st peak, the coordination number is increased almost exclusively to around 7, but considering both 1st shell and the intermediate pseudo shell, it becomes 10.6 with larger contributions from coordination numbers 3, 4, and 7. The range of second shell coordination numbers however is decreased from 26 to 37 in their earlier study to 16-28. As in their previous study, in this study also the Au + -water potential used for their MM simulation is not the widely used Lennard-Jones (LJ) potential, but a complicated combination of two-body and three-body terms. Both the above studies have dealt with Au + ions in water using QM/MM methods. Very recently, in an interesting study, ab-initio quantum mechanics (QM)-based Car-Parrinello molecular dynamics (CPMD) simulations have been used to study hydration behaviours of Au + ion and Au atom in bulk water [24]. It is interesting to note that they have not considered any negative counter ion of Au + ion. The Au + -Ow RDF (g(r)) obtained by this CPMD simulation [24] is strikingly different from that reported earlier [22,23]. In this case [24], a very sharp 1st peak followed by a "zero density" region between the first and the second peaks is observed. This region of zero density (or g(r)) is completely absent in the RDFs of both the QM/MM studies [22,23].
In the above-mentioned studies, computationally demanding quantum description of the force has been used either partially [22,23] or fully [24]. However, in order to study large systems like those involving membranes and other absorbing substances, quantum mechanical simulations of the likes of Car-Parrinello molecular dynamics (CPMD) and Born-Oppenheimer MD simulations are not practically viable due to their high computational costs and, therefore, use of classical molecular dynamics (MD) simulations [19][20][21] is one of the most suitable methods to explore the molecular-level properties, interaction mechanism, structure of a system of large size, and the dynamics at femtosecond resolution, over a longer period time, if accurate force fields for the system are known. In many previous studies [12,22,23], various potentials other than the Lennard-Jones (LJ) potential have been used. But most of the FF-based MD simulation methods as implemented in the available standard MD simulation program packages use pair potentials of the LJ type. It is important to mention that the classical LJ FF-based MD simulation has been used by Heinz et al. [15,16] to study gold atoms as face-centred cubic (FCC) metal. In these classical MD simulation studies, they have used 12-6 LJ and 9-6 LJ-type potentials to model interaction among neutral Au atoms in their solid state. The bulk properties such as density and surface tension of water near gold FCC metal have been predicted accurately. In another study, instead of Au + and Cl − ions in water, [AuCl 4 ] − complex ions in water have been studied [12] by using classical FF-based simulation, but the potential model used here is again a complicated one.
Since QM-based CPMD [24] and even the hybrid QM/MM simulations [22,23] are not viable for handling large lengthscales with appreciable timescales and results from the existing classical models are unable to produce results closer to the CPMD results, in the present investigation, we, therefore, intend to use a modified LJ potential for Au + -water interactions in classical MD simulations such that better results on the structural and dynamical properties of the aqueous solution of Au + ions can be obtained. One of our intensions is to use the LJ potential so that it can be amenable to any of the standard MD simulation program packages [25][26][27][28] so that large and complicated systems involving Au ions can be simulated easily. By using the modified LJ potentials, here, we obtain results closer to CPMD results as compared to the existing classical and QM/MM simulation results. It reproduces the peak position, height, and the appearance of the 'zero density' region between the 1st and 2nd peaks of the Au + -OW RDF in closer agreement with those obtained from the CPMD simulation [24]. We also characterize the first solvation shell in terms of geometry and orientation of water molecules around Au + ions and also study the residence time of the water molecules in the first solvation shell of the Au + ion.
In what follows, in Sect. 2, we describe Models and Methods, followed by Results and Discussion in Sect. 3. We analyse the simulation results of one Au + ion in water to obtain the geometrical structure and nature of the solvation shell and then we describe the residence time of water in the solvation shell of the Au + ion. A brief Concluding Remark will be offered in Sect. 4.

Lennard-Jones Potential Parameters
The gold ion is modelled as singly (positive) charged spherical ion interacting with each other and with OW of water through LJ potential [29,30], and the functional form of which can be written as where V LJ is the LJ interaction potential between two sites (denoted as i and j) separated by a distance of r ij . The LJ parameters ij and ij refer to usual energy and size parameters of the LJ potential. In the present work, we have used the units of kJ·mol −1 for ij and Å for ij . Many simulations with different sets of σ and ε parameters for the Au + ion and the SPC/E [31] water model have been performed to get the correct set of σ and ε so that same nature of the RDF of Camellone et. al. [24] can be reproduced. The details of the effects of variations of σ and ε parameters on the properties of ion-water RDF and coordination/hydration numbers have been displayed in Tables S1 to S6 in the Supporting information (SI). Apart from the SPC/E water model, other water models like Flexible-SPC and TIP3P have also been tried and the results from these water models have been found to produce similar ion-water RDFs as that obtained by the SPC/E water model (See Fig. S1 of the SI). In the present study, we have used the well-known SPC/E model [31], non-bonded potential of which is also of LJ form and this model is known to reproduce properties of bulk water at ambient temperature exceptionally well. For the interaction between a gold ion and water molecules, the cross interaction parameters for the LJ potential are obtained by using Lorentz-Berthelot (LB) mixing rule [30] i.e. arithmetic mean for ij and geometric mean for ij .

Computational Methods
All classical MD simulations were performed using the standard MD simulation package GROMACS, version 4.6.1 [25,[32][33][34]. A cubic simulation box of size 84.299 Å × 84.299 Å × 84.299 Å containing 1 Au + ion and 1 Cl − ion and 19998 SPC/E water molecules has been used. For better statistics, we have also simulated a system with 10 Au + and 10 Cl − ions in a box of 19980 water molecules. For comparison, we have also simulated similar systems in which non-bonded parameters for the Au + ion are obtained from the FF used by Heinz et al. [15] with a+ 1 charge on the Au atom.
In the NVT ensemble, the temperature of the system has been fixed at 300 K using V-rescale approach [35]. In case of simulations in NPT ensemble, the pressure of the system has been maintained at 1 bar using the Berendsen pressure coupling scheme [36][37][38]. The equations of motion have been integrated using velocity Verlet algorithm [19] with a time step of 2 fs. A real-space cutoff of 1.4 nm is used for calculating van der Waals interactions. The PME scheme [39,40] with a real-space cutoff of 1.4 nm is used to treat all Coulomb interactions. The bonds and angles of the water molecule are restrained to be fixed by using SETTLE algorithm [41]. Before starting a NPT run, a 2 ns simulation run in the NVT ensemble was performed. After that, a 2 ns NPT simulation run was used for equilibration, followed by a NPT production run of 10 ns.

Structural Quantities
In order to analyse the structure of the solvation shell of the Au + ion in water, we have calculated various RDFs e.g. g Au-OW (r) and g Au-HW (r) and adopted the nearest neighbours analyses [42][43][44][45][46][47][48]. Moreover, we have calculated various spatial and angular distributions of the nearest neighbours [45][46][47]. In order to analyse the orientation of a water molecule in the first or second solvation shell with respect to Au + -OW distance vector, we have calculated distribution of angles that various molecular vectors of water, namely, dipole moment (DM) vector, -OH bond vector and the vector ⊥ to the plane of the water molecule or plane ⊥ vector, make with the Au + -OW distance vector.

Residence Time Correlation Function
In order to assess the stability or subtleness of the hydration shell, average residence time of the water molecules in the first solvation shell is a very important parameter. Residence time correlation function (RTCF) portrays average time decay or relaxation of a water molecule initially present in the solvation shell or in any confined space. The RCTF, R(t), is then given by [48][49][50] where the function θ i (t 0 + t) is 1 if the ith water molecule is present in the first solvation shell of the ion at time t 0 and continues to be there till time (t 0 + t) and zero otherwise. The angular brackets with subscript 0 in the above equation signify averages over time origins t 0 and N Shell is the number of water molecules in the 1st solvation shell at time t 0 .

Lennard-Jones Parameters
Several simulations of Au + -water systems with varying ε and σ parameters for the Au + ions and SPC/E model for water have been carried out. Various features of the Au + -Ow RDF, g Au-Ow (r), as obtained from different choices of ε and σ for gold ion are compared with the same obtained by using the ab-initio CPMD simulation [24,51]. By comparing the data on peak positions, peak heights and coordination numbers as obtained from these simulations with the same obtained from ab-initio CPMD simulation [24,51], the most suitable non-bonded parameters for Au + ion in its aqueous solution have been found to be σ = 2.627 Å and ε = 0.1 kJ·mol −1 . Therefore, results presented in rest of this paper are based on this set of non-bonded parameters for Au + ions along with the well-known SPC/E model [31] for water with the cross parameters for Au + -Ow interaction obtained according to LB mixing rule [30].

Radial Distribution Function and Coordination Number
The RDF of oxygen and hydrogen atoms of water around Au + ion for our prescribed set of LJ parameters is shown in Fig. 1a. We have plotted the same obtained from our simulations along with those obtained using the modified (see Computational Methods section) FF of Heinz et al. [15,16] in Fig. 1b and the same obtained [51] from Ab-initio Car-Parrinello molecular dynamics (CPMD) simulations of Camellone et. al. [24] in Fig. 1c. It can be observed that the peak heights and peak positions of the Au-OW and Au-HW RDFs as obtained from the simulations using our newly proposed FF for the Au + ion are in close agreement with the work carried out by Camellone et. al. [24] using ab-initio CPMD calculations. Comparing the results (see Fig. 1a) obtained from our proposed FF with the results shown in Fig. 1b [15] and 1c [24], it is evident that our results (Fig. 1a) obtained by using our proposed FF parameters are closer in agreement with the ab-initio results (Fig. 1c) than the other results (Fig. 1b) obtained using another classical FF [15]. The most striking feature of the CPMD-derived Au + -Ow RDF is the appearance of the zero-density region between the first and the second RDF peaks. All the previous RDF results whether based on classical FF simulation [15] or QM/MM [22,23] simulations do not have this feature.
However, our new FF-based RDF (see Fig. 1a) shows a "zero density" region, in which g(r) is almost zero. Another important feature of the CPMD generated RDFs that the Au + -HW RDF starts from the tail end (minimum) of the 1st peak of Au + -Ow RDF is reproduced almost correctly in our work, but not in the earlier FF-based 15 and QM/MM-based [22,23] studies. The detailed features of the Au + -Ow RDF along with the coordination numbers of the 1st and the second hydration shells of Au + ion are given in Table 1. The 1st peak position and peak height of the Au-OW RDF as obtained from the newly proposed FF are 2.3 Å and 8.7 and the corresponding CPMD results are 2.1 Å and 8.6. The present results on peak position and height are in fairly Fig. 1 Au + -OW (red lines) and Au + -HW (blue lines) radial distribution functions obtained from a our present force field b force field taken from Ref. [15] and c CPMD simulation [51] (Color figure online) good agreement with the same from QM-based CPMD simulation [24]. In passing, it is worthwhile to mention that considering multiple Au + ions instead of one ion in water does not change the nature and details of the ion-water RDF (see Fig. S2 of the SI). However, as far as coordination number is concerned, all the FF-and QM/MM-based methods including the present FF-based one fails to reproduce the 1st shell coordination number of 2 as obtained in the CPMD simulation of Camellone et. al. [24] (see Table 1). It is important to emphasize that although QM methods, unlike the FF-based methods, take into account the electronic degrees of freedom, many QM-based ab-initio MD simulations cannot describe even the RDF of bulk water [52] correctly. In the absence of any experimental data on the RDF, of the same system, it therefore remains inconclusive whether QM-based MD simulation is superior to the FF-based MD simulations in describing the correct hydration structure of the Au + ion.
Apart from average coordination number, exact distribution of coordination number of the Au + ion is also an important characteristic, which tells us the proportion of ions with variable coordination in the solvation shell. Therefore, in Fig. 2, we have shown the distributions of first nearest neighbours P(N wat ) 1st shell and second nearest neighbours, P(N wat ) 2nd shell . For the 1st solvation shell, it is observed (see Fig. 2a) that the overall distribution is of non-Gaussian nature, and skewed around 6, indicating that in most of the cases, Au + ion is only 5-or 6-coordinated. To be specific, in almost 65% of the configurations, the ion is 6-coordinated, and in more than 30% configurations, it is 5 coordinated, whereas populations of all other (1-, 2-, 3-, 4-, 7-coordinated) species being almost zero. However, the distribution for the second shell neighbours is fairly symmetric with Gaussian nature, implying the presence 14 to 24 s shell neighbours (see Fig. 2b) with the peak at around 19. It is important to note that the corresponding results for classical simulation with the FF taken from Ref. 15 as shown in Fig. 2c are fairly symmetric with significant populations of 5-, 6-, 7-, 8-and 9-coordinated Au + ions, with a peak in the distribution at around 7. The peak of the 2nd shell distribution for the same FF (see Fig. 2d) is around 21, which is more than the same obtained from our proposed FF (cf. Figure 2b). Our result on second shell average coordination number of 19.2 (see Table 1) is closer to that (average coordination number 18) obtained from CPMD simulation (see Table 1).

Structure/Geometry of the First Solvation Shell of Au + Ion in Water
Having obtained average coordination number of 5.6 and knowing that mostly there are 5and 6-coordinated Au + ions for our newly designed Au + FF, we shall now attempt to establish the exact spatial arrangements of the neighbouring water molecules around the Au + ion in water. Various quantities like average neighbour distance and its fluctuation along with the distributions will give a rough idea about the solvation shell. Figure 3a shows the average distances of different neighbouring water molecules of the Au + ion along with its fluctuation (see inset). For this purpose, we have sorted the neighbours in terms of their distance from the Au + ion (identity of the water molecule may vary). It is interesting to observe that the average distances of the first four (distance wise) neighbours are around 2.3 Å (note that the first g(r) peak is at around 2.3 Å), but the same for the 6 th neighbour is around 3 Å. More importantly, the nature of the curve is sigmoidal, which indicates that there is sudden fluctuation in average distances for some of the neighbours and it also signifies that the first and the second solvation shells are distinct. Beyond the 6 th neighbour, average distances of the neighbours again vary within a small range of 3.8 Å to 4.8 Å and these molecules constitute the 2nd solvation shell (note that 2nd peak in g(r) appears at around 4.5). This is nicely captured by calculating the root mean square deviation, σ (fluctuation), of the average distances as shown in the inset of Fig. 3a. It is evident from this figure that the 6th neighbour has the maximum fluctuation (σ), but for all other neighbours, whether they are in the 1st and 2nd solvation shells, σ is very small. The large fluctuation (see inset of Fig. 3a) in the average distance of the 6 th neighbour clearly indicates that the sixth water molecule is at the boundary of the first and the second solvation shell and goes in and out of the 1st shell. A better picture of the relative position of the water neighbours with respect to Au + ion will be obtained if we calculate the distribution of distances of different neighbours. In Fig. 3b, we have, therefore, shown the distribution of n th water neighbour distance from the Au + ion. As clearly seen in this plot (Fig. 3b), the first five neighbours have the unimodal distributions with the spread of 2-3 Å. The unimodal distribution of distances of the 7th or 8th neighbour with the peak of the distribution at around 4 Å is a confirmatory signature that these neighbours belong entirely to the 2nd solvation shell of the Au + ion. It is the 6th neighbour, which has a bimodal distribution with the two peaks around 2.5 and 3.8 Å, indicating that in a significant number of configurations, it goes out of the 1st solvation shell and enters the 2nd solvation shell. The fact that sixth neighbour is at the boundary of the 1st solvation shell and it toggles between the 1st and the second solvation shells is also evident from the appearance of a maximum in the plot of fluctuations in average distances of different neighbours (see inset of Fig. 3a). It is interesting to observe that the results from our simulation with the FF taken from Ref. 15 are presented in Fig. 3c and d, and the average distance of the neighbours is increasing continuously without any sudden jump, which is observed in Fig. 3a. Further, by plotting the distributions of distances in Fig. 3d, unimodal distributions for all the eight nearest neighbours with continual shifts in peak positions have been found. It is important to note that such features as the sigmoidal nature (cf. Figure 3a) of average distances of the neighbours and bimodal nature (cf. Figure 3b) of the distribution of distances for some of the neighbours are not observed in case of results obtained from the FF of Ref. [15] (see Fig. 3c and d). It emphasises that such intricate nature of the solvation shell depends on the critical choice of the force-field parameters.

Ow-Ow (inter-neighbour) Distance Distributions for First Eight Water Neighbours
In order to further confirm the octahedral nature of the 1st solvation shell, we have calculated distributions of distances between any two water neighbours. In this calculation, we have considered first 8 neighbours (2 more than 6 neighbours forming octahedral structure) and calculated their distances with 1st six first shell water neighbours. These distributions are shown in Fig. 4. According to the octahedral arrangement of 6 water neighbours around the central Au + ion, it is expected that if the Au + -OW distance (assuming Au + -OW distance be same for all first six neighbours) be a then we should expect two peaks, namely, one at a √ 2 and another at 2a , in the distributions of intermolecular distance between any two of the first six water neighbours. From the 1st peak position of the Au + -OW RDF, a can be estimated for our case to be 2.3 Å. Therefore, we can expect peaks at around 3.25 Å (= a √ 2 ) and another at 4.6 Å (= 2a ), and indeed we observe these two peaks in the distributions of first six neighbours as shown in Fig. 4. The distributions calculated for 7th and 8th neighbours show two peaks at around 2.7 Å and 5.5 Å. These are the distances of the second shell molecule with the 1st shell neighbours. Assuming the 2nd shell neighbour is approaching the 1st shell from a direction which bisects the OW-Au + -OW angle formed by two 1st shell water (OW) neighbours, and if the closest approach between the 1st and the 2nd shell water be around 2.7 Å, the second peak should arise at around 5.6 Å. We indeed find two such peaks at 2.7 Å and 5.5 Å in Fig. 4.
Having almost identified the octahedral nature of the first solvation shell of Au + ion from different distance distributions of nearest neighbours, more confirmatory evidence of the octahedral arrangement of the first six water neighbours can be obtained from the distribution of angles made by any two water neighbours extended at the Au + site.

Distributions of OW-Au-OW Angles for the 1st Solvation Shell Water Molecules
The distribution of angles extended to Au + ion by any two nearest water neighbours will give a clear idea about structure or geometrical shape of the first shell neighbours. We have identified first 8 neighbours and for each neighbour, we have calculated the above angle and its distribution, which is shown in Fig. 5. While calculating the angle for the 7th or 8th neighbour, the other water molecule chosen is any of 1st to 6th neighbours. It is interesting to observe that there are two prominent peaks in the Ow-Au-Ow distribution of first five neighbours at 90° (cos θ = 0.0) and 180° (cos θ = −1.0). The presence of these two peaks signifies that the first solvation shell structure is either trigonal bi-pyramidal (5 neighbours) or octahedral (6 neighbours). It is also seen (see Fig. 5) that the 6th neighbour Fig. 4 Distribution of distance between nth water neighbour and any other water neighbour from the first six nearest neighbours. For 7th and 8th neighbours, the distance is between that (7th or 8th) neighbour and any of the first six nearest neighbours 1 3 is also maintaining these positions mostly. In order to visualize the neighbour arrangement around the Au + ion, we have also analysed the trajectory and two types of solvation shells are observed. In some cases, a central Au + ion has five neighbours and in some others there are six neighbours. The snapshots of these two representative solvation shells with 5 and 6 neighbours are shown in Fig. 6a and b, respectively. From Fig. 6, it is evident that the Au + ion with 5 water neighbours forms a trigonal bi-pyramidal (slightly distorted) structure, whereas that with six neighbours forms an octahedral arrangement. Square pyramidal structure is unlikely as in that case one side of the Au ion will have no neighbour, which is an impossibility considering the spherical nature of the Au ion. In fact, by searching the entire trajectory, no such structure has been found.
The distributions for the 7th and the 8th neighbours are completely different from those of the first six neighbours. The distribution has peaks at around 45° (cos θ = 0.7) and 138° (cos θ = − 0.75). These molecules are in the 2nd solvation shell and their positions are complimentary to the neighbours in the first solvation shell to maximize the hydrogen bonding with the 1st shell neighbours. It also justifies the inter-neighbour distance distribution peaks for 7th and 8th neighbours (cf. Figure 4). Thus, from these analyses, we are in a position to predict the trigonal bi-pyramidal/octahedral geometry of the first solvation shell. The same distributions obtained from our simulations with the FF taken from Ref. For 7th and 8th neighbours, the angle is extended by that (7th or 8th) neighbour and any of the six nearest neighbours. Panel a shows the results from our FF and b the same calculated using FF taken from Ref. [15] Fig. 6 Snapshots of the solvation shells of the Au + ion (yellow colour) with a five and b six neighbours. The line connecting the Au + ion and the oxygen of water molecule is not showing the bonds. It is to guide the eye [15] are shown in Fig. 5b. In this case, we have seen (Fig. 5b) peaks corresponding to 72° and 140° confirming a pentagonal arrangement of the water neighbours around the Au + ion.

Orientation of Water Molecules Around Au Ions
Having established the almost octahedral structure/geometry of the 1st solvation shell, it is now important to know the orientation of the water molecules in the first solvation shell. A water molecule has three orientational vectors, namely, the dipole moment vector (DM), -OH bond vector and the vector perpendicular (⊥) to the molecular plane of water. In the following sub-section, we shall describe the distributions of different molecular orientational vectors of water with respect to the Au + -OW distance vector.

Angular Distributions of the Water Molecules in the First Solvation Shell
The schematic representations of the angles θ, ψ and φ made by molecular orientational vectors DM, -OH and plane ⊥ vectors, respectively, with the Au + -OW distance vector of water are shown in Fig. 7. The distributions of the corresponding angles (θ or ψ or φ) for the water molecules in the 1st solvation shell of the Au + ion are shown in Fig. 8. Here, the angle is the average of all angles made by all the 1st shell water neighbours. As evident from Fig. 8, the peak corresponding to the DM vector is at θ = 0° (see green coloured line in Fig. 8a), indicating the alignment of DM vector with the Au + -OW distance vector. Therefore, according to this plot, the H atoms of water are extended outward from the Au + ion. The DM vector is the angle bisector of molecular angle HW-OW-HW of water starting from OW of water. Since DM vector is at an angle θ = 0° with the Au-OW distance vector, it is therefore expected that the angle made by the -OH bond vector will have an angle of around half of the bond angle i.e. ~ 52° with the Au-OW distance vector, which, in turn, is aligned with the DM vector of the water. Therefore, the -OH bond vector should have a peak in the angular distribution around ψ = 52° or cos ψ = 0.6 and we observe a peak around cos ψ = 0.55 region in the distribution of angle ψ (see the red coloured line in Fig. 8a). With these two orientations of the water molecule fixed, a water molecule can rotate around the DM vector axis. Therefore, the ⊥ vector (which is perpendicular to molecular plane of water) can have any angle with the Au-OW distance vector. So, the distribution of angle made by the ⊥ vector with the Au-OW vector (see blue line in Fig. 8a) is flat with no preferred peaks in it. The same angular distributions obtained from  Fig. 8b. The results are almost the same as in Fig. 8a. Such orientation profiles have also been calculated [53,54] in case of aqueous solution of uranyl ions.

Neighbour-Wise Angular Distributions of the Water Molecules in the First Solvation Shell
In order to get neighbour-wise orientational distributions, we have calculated these distributions for all 1st shell water neighbours of the Au + ion and plotted them in Fig. 9. In all the distributions, up to 5th neighbours, we find almost the same orientation as the average one (displayed in Fig. 8). The sixth neighbour, as it is toggling between 1st and 2nd solvation shells, its time averaged orientation is slightly different from the others (see Fig. 9).

Residence Time Correlation Function of Water around Au + Ions
The mean residence time of a water molecule in the first solvation shell of Au + ion can be obtained from the residence time correlation function as calculated according to Eq. (2). The RCTF of the first solvation shell (R(t)) water molecules around the Au + ion is displayed in Fig. 10. The decay is purely single exponential indicating a uniform motion for all the water neighbours. The average or mean residence time (τ) of the first shell water is obtained by fitting the R(t) data by a single exponential of the type: where R(t) are the fitted RTCF data (see the open circle points in Fig. 10). The value of mean residence time (τ) calculated in this way is 17.36 ps. Rode and co-workers in their QM/MM study [22] in which only the first solvation shell was treated quantum mechanically, have found a very labile 1st solvation shell with τ = 7.5 ps. Because of the very small solvation time obtained in that study, they have thought it reasonable to include the 2nd shell also in the QM region. Accordingly, in their follow-up study, in which both first and second solvation shells have been treated quantum mechanically, they have found the mean residence time to be τ = 10.47 ps [23]. Our value τ = 17.36 is more than that such that the solvation shell is found to be reasonably stable, with a little exchange of water molecules between the first and second solvation shells. This result is commensurate with the 'zero density' region obtained between the 1st and the 2nd peaks of the Au-Ow g(r) (cf. Figure 1a). Residence times for other univalent ions as obtained from the classical FF-based simulations have been found [55] to be of the same range.

Conclusions
In summary, we have used a modified LJ-based force field for gold ions in SPC/E water. Using this modified FF, we have calculated various features of the Au + solvation shells. The Au + ion in most of the configurations is either six or five coordinated. We have used nearest neighbour analyses to find the geometry of the solvation shell. By analysing the distributions of Au + -OW and OW-OW distributions of the first shell water molecules, the geometry is roughly identified as octahedral (for hexa-coordinated) or trigonal  bi-pyramidal (for penta-coordinated). Further analysis of the distribution of OW-Au + -OW angles reveals the presence of 90° and 180° angles corroborating the octahedral and/or trigonal bi-pyramidal geometry of the solvation shell of the Au + ions. In order to have an idea about orientations of the solvation shell water molecules with respect to Au + -OW distance vector, we have calculated distributions of angles made by any of the three molecular orientational vectors, namely, dipole moment vector, -OH bond vector or the vector ⊥ to the plane of the water molecule, with the Au + -OW distance vector. The peak in the DM distribution appears at 0°, whereas that in the -OH distribution appears at around 50-52°, indicating that the water molecules are arranged with their H atoms away from the Au + ion. From the residence time correlation function, the average residence time of a water molecule in the first solvation shell of the Au + ion is found to be around 17 ps. The present FF can reproduce the peak height and zero-density region between the first and the second peaks of the Au + -OW RDF as observed in the CPMD simulations. It is also to emphasize that the present results are not only better than those obtained from other available classical force-fields-based simulations, but better than OM/MM results as well. The LJ-based force field for Au + ion presented in the present study will be amenable to any standard force fields included in the available MD simulation packages [25][26][27][28] and enable us to simulate large and complex systems for appreciable amount of time. Uses of this newly developed FF in investigating adsorption of the Au + ions in various porous materials at varying thermodynamic conditions are in progress.