Finite temperature dynamics of anionic water trimers
FFinite temperature dynamics of anionic water trimers
J. M. Finn a , F. Baletto a, ∗ a Physics department, King’s College London, Strand, London, UK, WC2R 2LS
Abstract
Utilising Car-Parrinello molecular dynamics simulations, the finite temperature behaviour of water trimers,(H O) , in both neutral and anionic frameworks, has been investigated. A significant structural change inthe anionic structure has been observed at temperatures above 100 K where a chain geometry has formedand stabilised entropically. On the other hand, neutral trimers have remained in their ring structure, aspredicted theoretically at low temperatures, for long time periods. Keywords:
Water trimer, Excess electron, Finite temperature
1. Introduction
Water clusters have a vast role within various areas of science, ranging from radiation chemistry [1] andmolecular biological processes [2, 3] to the suggestion they play a non negligible role in atmospheric chemistry[4, 5]. The hydrated electron, usually denoted as e − aq , in water systems represents an intriguing and importantexample of electron solvation, since its discovery in bulk water, where the hydrogen-network (H-network)of six water molecules form a hydration shell with an octahedral symmetry [6, 7]. Recent experimental andtheoretical works have demonstrated that a dynamical rearrangement of the surface free OH-groups occursboth at the ice-vacuum and water/air interfaces to localise an excess electron (EE) [8, 9, 10, 11, 12, 13].In the last two decades, thanks to the advances in cluster science, significant research has been undertakeninto the description of electron hydration in water nanoparticles and how to extrapolate this information tothe macroscopic level. Specific challenges include the minimum number of water molecules and geometricstructure required to bind an excess electron, which are found experimentally through electron injectiontime of flight spectroscopy and electron photodetachment spectra [14, 15]. It has been shown that thevertical detachment energy (VDE), which is the vertical binding energy of the excess electron evaluated atthe anionic cluster morphology, is related to the localisation/delocalisation of the extra electron, althoughit does not provide any information regarding the “position” of the excess electron with respect to thehydrogen network, i.e. surface or cavity bound. Besides several electronic structure calculations devoted tofind the best geometrical isomer for anionic water clusters, to estimate their VDE and to determine whetherthe extra electron is surface or internally bound [12, 16, 17, 18], the current understanding of how subtlethermal fluctuations of the aqueous solution influences the diffusive properties of the EE is still poor, dueto the variety of structural motifs that a water cluster has access to during its dynamics. Recently, through ab-initio Born-Oppenheimer molecular dynamics simulations it has been predicted that even on mediumsize water clusters, the EE is preferentially surface solvated [19]. Moreover, free jet ion source experimentshave shown that the VDE depends on the different carrier gas used, and thus on the cluster temperature,allowing the formation of different morphologies [20, 21]. It is of the utmost importance to have an atomisticdescription of the finite temperature dynamics of the H-network in order to give insights into the electrontransfer processes [22, 18].In this manuscript for the first time, the linearisation process of the anionic water trimer has beenobserved through finite temperature Car-Parrinello molecular dynamics. Although a linearisation mechanism ∗ Corresponding author: [email protected]
Preprint submitted to Elsevier October 25, 2018 a r X i v : . [ phy s i c s . a t m - c l u s ] F e b as previously been explored both experimentally [23], with an estimation of the attachment energy of 142 ±
2. Computational Methods
To study the finite temperature behaviour of the water trimer, in both the neutral and charged frame-works, at atmospherically relevant temperatures Car-Parrinello molecular dynamics simulations [30] havebeen performed, using the open source Quantum-Espresso package [31]. The electrons were attributed afictitious mass of 400 a.u. and the Brillouin zone sampling was restricted to the Γ-point only. A timestepof 1 fs (4 a.u.) was used to integrate the equations of motion. The valence electron-nuclei interactions weredescribed by Troullier-Martins pseudopotentials [32]. Becke-Lee-Yang-Parr (BLYP) exchange-correlationfunctional [33] has been considered to describe the valence electron-nuclei interactions with a plane waveenergy cut-off of 100 Ryd, to enable the system to be computationally tractable and the energy convergencethreshold is less than 0.1 mHa. Periodic boundary conditions have been applied to a cubic box of length40 bohr, and a multiple vibrational Nos´e-Hoover thermostat [34], characterised by two frequencies of 8 THzand 20 THz to thermalise both the slow oxygen-oxygen and the fast oxygen-hydrogen motion, respectively,has been used. The system has been heated gently to around 50-80K and allowed to equilibrate before beingheated again to 120K. The whole thermalisation time was at least 5 ps and the reported data have beenaccumulated over 15 ps for both the neutral and charged water trimers, with over 100 ps obtained for theneutral trajectory and 25 ps for the higher temperature dynamics.The dynamical analysis of the charged system has been performed within the Mauri self-interactioncorrection scheme (MSIC) [35], which has been shown to give good results for the description of watersystems [8, 9, 36]. This framework works in the local spin density approximation where the paired spinup and spin down wave-functions are forced to be the same and the energies associated with the Hartreepotential, E H , and the exchange and correlation functional, E xc , are corrected such that; E M,H = E H [ n ↑ + n ↓ ] − (cid:15)E H [ n ↑ − n ↓ ] E M,xc = (1 − α ) E xc [ n ↑ , n ↓ ] + αE xc [ n ↓ , n ↓ ] , (1)where n ↑ = n ↑ ( r ) and n ↓ = n ↓ ( r ) are the charge densities for spin up and spin down, respectively. We haveused a modified version of the original MSIC scheme whereby the self-interaction is only partially screened,imposing that the two parameters α and (cid:15) , introduced by Sprik and coworkers [37], are both less than theunit. The choice of these two parameters follows the requirement to reproduce the electron attachmentenergy of a water dimer of 26 −
51 meV [38]. Using α = 0 .
85 and (cid:15) = 1 .
00, a corrected value of 49 meV hasbeen calculated against an enormous binding energy of 911 meV using an uncorrected local spin density.Due to the imposed periodic boundary condition, a Madelung correction [39, 40] has to be added in order2o remove the spurious Coulomb interaction among periodic images. This contribution is proportional to(1 − (cid:15) ) because the interaction of the singular occupied molecular orbital with itself has been already partiallyremoved when working in the Mauri scheme.
3. Results
The water trimer presents different isomers which can be divided in two large classes: cyclic and linear.A structure is labelled as ring or cyclic when three hydrogen bonds with length between 1.2 and 2.7 ˚A areformed. To be classified as linear there are two requirements that have to be satisfied simultaneously: aninter-oxygen distance R OO larger than 3.8 ˚A; and a obtuse internal angle, Θ OOO . Among ring clusters afurther characterization has been carried out throughout the description of the geometrical position of thedangling OH moieties, accordingly to their direction perpendicular to the oxygen plane, where up ( u ), down( d ) are arbitrarily defined as opposite to each other, and planar ( p ) means that the OH lies on the oxygenplane. The planar position is where the free OH group lies in a position ± .
05 ˚A from the oxygen plane.Any groups greater or smaller than this are labelled as up and down, respectively.This classification is the basis of the proposed Isomer Population Analysis (IPA) which allows the de-termination of the cluster morphology, including linear shapes, at each frame in the dynamics in order tocalculate the frequency of each visitation.For the neutral case, the most stable structure has been calculated to exist as a cyclical hydrogen network,where two of the dangling OH groups point in the same direction with respect to the oxygen plane, andthe other in the opposite direction ( uud ), in agreement with other high level calculations [41, 42]. The nextfavourable structure is the uuu where all the free OH point in the same direction, and the final ring consideredis the upd which has been found to be the transition state between two degenerate uud isomers. Each ofthe six equivalent uud isomers are, indeed, easily accessible through an external hydrogen flip through theoxygen plane. The activation energy barrier between two neighbouring uud neutral water trimer is of theorder of 38 meV, calculated by means of nudged elastic band [43]. On the other hand, the transition fromthe uud to the uuu structure presents an activation energy barrier of 68 meV.A neutral cluster, starting from a uud isomer, has been allowed to evolve at 120 K for 100 ps and theanalysis of its dynamics is reported in blue boxes in Figure 1. No structural changes from the cyclical isomershave been observed, the uud remains stable for around 50% of the total simulation time, with a only a shorttime interval (10 %) spent in the less energetically favourable uuu geometry. A frequent hydrogen motionhas been detected and the transition between two equivalent uud rings has been observed in agreement withthe nudged elastic band calculations of hydrogen flipping. It is worth stressing that at no point in time isthe linear structure observed, as the energy loss due to the breaking of a hydrogen bond is much too high( ∼ uuu isomer and has been allowed to evolve for 10 ps,where like the neutral dynamics, no structural transition have been detected from the ring series of isomerswith the uud shape found more frequently than the uuu , although the latter has a higher total dipolemoment. Over the total time period the internal angles are not observed to change further than 60 ◦ ± ◦ ,and the ionic temperature remains relatively stable. The IPA shows that for over 80% of the simulationtime the uud isomer is found, with a further 9% in the uuu , as labelled by red boxes in Figure 1.The second anionic molecular dynamics has been performed at a higher temperature, again starting fromthe uuu isomer, and heated using a series of thermostats to around 120 K and then left free to evolve withinthe Car-Parrinello scheme. The total thermalisation time including the thermostats was around 8 ps, whichwere discounted from the analysis. Within 2 ps after the equilibration one of the hydrogen bonds brokecausing an increase in the maximum internal oxygen angle, indicating the ‘opening’ into the linear isomer.The chain isomer is observed to last for around 15 ps, when the ring series of isomers are reformed andremain stable for an additional 2 ps. IPA analysis of the trajectory shows that the linear isomer is dominantfor around 62% of the simulation time, with a secondary isomer of the uud structure, as depicted by greenboxes in Figure 1. 3wo distinct mechanisms occurred during the total dynamics which allowed the water trimer to tem-porarily remain stabilised in this maximal dipole form. The first, highlighted by the shadow region in Figure2 labelled with Stage I , is the breaking of the hydrogen bond, thanks to a rotation of a free OH groupwhich causes torsional stress on the bond itself. This allows the water molecule, previously donating itshydrogen, to rotate away from the acceptor molecule, and freeing both OH groups. However, the cluster isstill not stable and large oscillations of the internal angle cause the isomer to attempt to recombine into thering structure on the very short time scale, as reported in the bottom panel of Figure 2. The second mecha-nism, grey region on the bottom panel of Figure 2 titled
Stage II , describes a bifurcation flip, whereby thewater molecule is seen to rotate through a large angle perpendicular to the oxygen plane, initially breakingthe hydrogen bond being donated to the middle molecule in the linear chain. As the rotation changes thehydrogen network and the bond breaks, the second hydrogen moves into the same position in order to formagain the chain structure, thus the free OH groups are swapped.Static analysis of few structures, taken from the dynamical trajectory, gives some indication of the drivingforce behind the linearization process observed at the higher temperature. These structures are depictedin Figure 3, the cyclical uud (R ) and uuu (R ) motifs in the top row, and the linear shape where L was constructed by hand with an internal angle of 180 ◦ and L and L are configurations from the highertemperature dynamics at the end of STAGE I and STAGE II, respectively. Using the partial MSIC correctionof α = 0 .
85 and (cid:15) = 1 different snapshots have been analysed and VAE calculated (table 1). Correspondingcalculations using the B3LYP xc-functional confirm the results that the uud is the least bound structure,and L is the strongest. All of these results are comparable with literature values, using both MP2 and CIapproaches, with strengths of -60 meV and -243 meV respectively [44, 45].The spin density difference between the spin-up and spin-down charge densities population gives anindication of the position of the additional charge: isosurfaces have been plotted in Figure 3 depicting 70%of the total charge. In agreement with the VAE results both the linear structures taken from the dynamicsshow spatial localisation of the additional charge, at the double acceptor end of the chain. Moreover, it wasobserved that not only was the localisation and formation of the anionic structure due to the increase indipole moment of the total cluster but also due to the directionality of the individual dipole moments ofeach water molecule. The dipole moment of the uuu isomer, which is above the critical amount required tobind the additional charge, however, was crucially not found to attach the additional charge, but a structurewhich had a co-operative dipole moment formation such as L which has a smaller dipole does attach, withonly a small difference in the attachment energy. Additionally, once the second stage of the dynamics hadoccurred, and the final linear isomer was formed (L ) the additional charge became bound to the structurewith an attachment energy of -205.5 meV, and spatially was observed to form a bean shaped cloud near tothe double acceptor at the end of the chain. Thus this co-operative effect of the individual dipole momentsof each water monomer in the cluster is found to have a strong impact on the formation and localisation ofthe additional charge.To support the idea that the dynamics are driven through a kinetic attachment the vibrational densityof states (VDOS) has been calculated V DOS ( ω ) = (cid:18) √ π (cid:90) ∞−∞ e − iωt Z ( t )d t (cid:19) (2)where Z ( t ) is the velocity autocorrelation function Z ( t ) = (cid:104) v (0) · v ( t ) (cid:105)(cid:104) v (0) · v (0) (cid:105) . (3)Using the VDOS the entropic contribution for each of the finite temperature trajectories can be calculated[46, 47]: T ∆ S = 3 k B (cid:90) ∞ V DOS ( ω ) f ( ω ) dω, (4)4eing f ( ω ) = (cid:104) ¯ hω k B T coth( ¯ hω k B T ) − log(2 sinh( ¯ hω k B T )) (cid:105) where and k B is the Boltzmann constant and T is thetemperature. The change in entropy between the uud and linear structures within the anionic framework isfound to be -0.025 eV/K, showing that the linear form is entropically stable at higher temperatures, whichconfirms that there is a kinetically driven attachment of the additional charge.The VDOS for each trajectory has been plotted in Figure 4 where finite temperature dynamics in theneutral framework agree well with available calculations by Kang et al. [48], with each of the three bandsbeing expressed. The first band, 0-1000 cm − shows more peaks in the anionic dynamics correspondingto vibrations in the cluster during the opening and subsequent stabilisation of the chain isomer. Theattachment of the excess electron is suggested by the distinctive peak at 250 cm − seen during both anionicdynamics. The formation of a double-acceptor single molecule can be attributed to the different shape ofthe vibrational bend peak around 1500 cm − , in agreement with experimental observation [49]. However,for both anionic trajectories the third band, between 3000-4000 cm − corresponding to intra-molecular OHstretching, is under expressed.
4. Discussion
Car-Parrinello molecular dynamics simulations of both neutral and charged water trimers have been per-formed at temperatures relevant to atmospheric processes. Our simulations have shown a direct correlationbetween the directionality of the individual dipole moments of each water monomer and attachment andspatial localisation of the excess electron. Furthermore, it has been observed that when the anionic systemwas at a temperature above 100 K the linear formation of the cluster is dominant for around 62% of thetotal simulation time, with a temporal stability of around 15 ps.The fast transition from the ring to linear isomers at 120 K indicates that not only is the energeticstability important during the formation of a water cluster anion, but also that the total dipole, enablingthe binding and localisation of the additional charge plays an important role.The dipole decomposition has allowed the characterisation of the linear isomers into two types, the firstrelates to the period of time prior to the bifurcation flip, where the cluster attempts to form co-operativedipole moments in order to maximise the dipole along a vector towards the charge cloud. The second typeoccurs subsequent to the flip, where the dipole moments are able to align towards the additional chargeand therefore, localise the excess electron, co-operatively binding through multiple dipole moments pointingtowards the charge. These two effects appear to be behind the driving force of the dynamics.Thus, the co-operative dipole binding could have significant implications for larger water clusters, andcould potentially be behind the general method for binding an additional charge. Secondly an entropiceffect has been observed allowing the formation of the linear anion, and the attachment and localisationof the additional charge when the temperature is in excess of 100K. This implies that modelling the finitetemperature behaviour of the attachment of excess electrons is important in order to understand the bindingmotifs, and have greater applicability to the atmospheric anionic formation.
5. Acknowledgments
The authors are grateful to the UK research council EPSRC for its financial support. JF appreciate theHPC2-Europa and Cineca facilities for their financial and computational support. Authors would like tothank the King’s College London system manager, Dr. A. Comisso, for the management of local computa-tional facilities.
6. ReferencesReferences [1] B. C. Garrett et al.,
Chem. Rev. (2005) 355.[2] F. Garczarek, K. Gerwert,
Nature (2006) 109.[3] I. Michalarias, I. Beta, R. Ford, S. Ruffle, J. Li,
Appl Phys A-Mater (2002) S1242.
4] J. Headrick, V. Vaida,
Phys Chem Earth Pt C (2001) 479.[5] Q. B. Lu Phys. Rev. Lett. (2009) 118501.[6] E. J. Hart, J. W. Boag,
J. Am. Chem. Soc. (1962) 4090.[7] S. Schlick, P. Narayana, L. Kevan, J Chem Phys (1976) 3153.[8] F. Baletto, C. Cavazzoni and S. Scandolo, Phys. Rev. Lett. (2005) 176801.[9] U. Bovensipene et al., J. Phys. Chem. C (2009) 979.[10] A. Madarasz, P. J. Rossky, L. Turi,
J Chem Phys (2007) 234707.[11] D. M. Sagar, C. D. Bain, and J. R. R. Verlet,
J. Am. Chem. Soc. (2010) 6971 .[12] J. M. Herbert and L. D. Jacobson,
Int. Rev. Phys. Chem. (2012) 1.[13] B. Abel, U. Buck, A. L. Sobolewski, W. Domcke, Phys. Chem. Chem. Phys. (2012) 22.[14] J. Verlet, A. Kammrath, G. Griffin, D. Neumark, J. Chem. Phys. (2005) 231102.[15] O. T. Ehrler, D. M. Neumark,
Acc. Chem. Res. (2009) 769.[16] H. Tachikawa, J. Chem. Phys. (2006) 144307.[17] K. Yagi, Y. Okano, T. Sato, Y. Kawashima, T. Tsuneda, K. Hi- rao,
J. Phys. Chem. A (2008) 9845.[18] O. Marsalek and F. Uhlig and J. Vandevondele and P. Jungwirth,
Acc. Chem. Res. (2012) 23.[19] T. Frigato, J. VandeVondele, B. Schmidt, C. Schuette, P. Jungwirth, J. Phys. Chem. A (2008) 6125.[20] R. M. Young and M. A. Yandell and S. B. King and D. M. Neumark,
J. Chem Phys. (2012) 094304.[21] R. N. Barnett, R. Giniger, O. Cheshnovsk, and U. Landman,
J. Phys. Chem. A , (2011) 7378.[22] K. A. Tay and A. Boutin, J. Phys. Chem. B (2009) 11943.[23] L. K.Liu, J. Loeser, M. Elrod, B.Host,
J Am Chem Soc , (1994) 3507.[24] S. T. Arnold, PhD Thesis, Photoelectron spectroscopy of solvated electron and solvated anion clusters (1993), ChemicalDepartment, John Hopkins University.[25] A. D. van der Avoird, and K. Szalewicz,
J. Chem. Phys. (2008) 014302.[26] T. Taketsugu and D.J. Wales,
Mol. Phys. (2002) 2793.[27] P. Ayotte, et al.,
J. Chem. Phys. , (1999) 6268.[28] N.I. Hammer, and J.R. Roscioli, and M.A. Johnson, MA, J. Phys. Chem. A (2005) 7896.[29] F.N. Keutsch and R.J. Saykally,
PNAS (2001) 10533.[30] R. Car and M. Parrinello, Phys. Rev. Lett. (1985) 2471[31] P. Giannozzi et al., J. Phys-Condens. Mat. (2009) 395502.[32] Troullier, N and Martins, JL, Phys. Rev.
B43 (1991) 1993.[33] A.D. Becke,
Phys. Rev.
A38 (1988) 3098.[34] S. Nose,
J. Chem. Phys. (1984) 511.[35] M. dAvezac, M. Calandra, F. Mauri, Phys Rev
B 71 (2005) 205210.[36] S.Bhattacharya, J. Finn, V.Diep, F. Baletto, S. Scandolo,
Phys Chem Chem Phys (2010) 13034[37] J. VandeVondele, M. Sprik, Phys Chem Chem Phys (2005) 1363.[38] Y. Bouteiller, C. Desfrancois, J. Scherman, Z. Latajka, B. Silvi, J Chem Phys (1998) 7967.[39] G. Makov, M. Payne,
Phys Rev
B 51 (1995) 4014.[40] I. Dabo, B. Kozinsky, N. E. Singh-Miller, N. Marzari,
Phys Rev
B 77 (2008) 115139.[41] F. Keutsch, J. Cruzan, R. Saykally,
Chem Rev (2003) 2533.[42] T. Salmi, H. G. Kjaergaard, L. Halonen,
J Phys Chem A (2009) 9124.[43] G. Henkelman, H. Jonsson, J Chem Phys 113 (2000) 9978.[44] S. Iwata, F. Chen, J Elect Spectro 142 (2005) 277.[45] T. Tsurusawa, S. Iwata,
Chem Phys Lett (1998) 553.[46] U. F. T. Ndongmouo, M. S. Lee, R. Rousseau, F. Baletto, S. Scandolo,
J Phys Chem A (2007) 12810. 014302.[47] R. de Sousa, H. W. Leite Alves,
Braz J Phys (2006) 501.[48] D. Kang, J. Dai, Y. Hou, J. Yuan, J Chem Phys (2010)[49] J.Roscioli, N. Hammer, M. Johnson,
J. Phys. Chem. A (2006) 7517. VAE [meV]Structure µ [D] SIC-CP[0.85,1] B3LYPRing motifuud (R ) 1.24 -63.1 66.2uuu (R ) 4.05 -86.5 -0.04upd 0.98 -62.0 98.6Linear motifL Table 1: Vertical attachment energy ( E V AE in meV) and total dipole moments ( µ in Debye) for the considered structuresreported in Figure 3. Results are compared with the B3LYP functional. A negative value in the VAE refers to a bound EE.Experimental reference value is a bound state at 142 ± !" uud uuu upd ppd pdd upp ppp lin Isomer0 . . . . . T i m e [ % ] . . . . . Figure 1: Comparison of the IPA for both the neutral (blue) and anionic 80K (red) and 120K (green) finite temperature dynamicsof the water trimer cluster, clearly showing that for the anionic cluster at 120K the linear isomer is chosen preferentially, whereasfor the neutral and anionic 80K dynamics the ring isomers are the only structures found. The time is expressed as a percentageof the total simulation time and is unit-less. igure 2: Evolution of the internal oxygen angles within the water trimer at temperatures around 100 K. Top panel shows thetemperature evolution (K) and bottom panel shows the two internal oxygen angles ( ◦ ). Shaded regions indicate the two stagesof the dynamics, as described in the text.Figure 3: Isosurfaces depicting 70% of the total charge of the excess electron in differing water trimer isomers. In additiongreen arrows show the individual dipole moments for each water molecule: notice the uncooperative dipole moments in the uuu structure, and the cooperative binding allowing for the localisation in the linear isomers. N o r m a li s e dd e n s i t y o f s t a t e s [ a r b ] CHARGED (120K)CHARGED (80K)NEUTRAL v [cm − ] Figure 4: Vibrational density of states for the three water trimer dynamics, neutral and anionic. Frequency is in cm − andnormalised density of states is in arbitrary units. Neutral dynamics at 120 K.andnormalised density of states is in arbitrary units. Neutral dynamics at 120 K.