Transcription
Articlepubs.acs.org/JPCBAb Initio Molecular Dynamics Simulations of an Excess Proton in aTriethylene Glycol Water Solution: Solvation Structure, Mechanism,and KineticsMarshall T. McDonnell,† Haixuan Xu,‡ and David J. Keffer*,‡†Department of Chemical and Biomolecular Engineering, and ‡Department of Materials Science and Engineering, University ofTennessee, Knoxville, Tennessee 37996, United StatesS Supporting Information*ABSTRACT: We investigate the solvation shell structures, thedistribution of protonic defects, mechanistic details, kinetics, anddynamics of proton transfer for an excess proton in bulk waterand for an excess proton in an aqueous solution of triethyleneglycol (TEG) via Car Parrinello molecular dynamics simulations. The PW91, PBE, and PBE with the Tkatchenko Scheffler (TS) density-dependent dispersion functionals wereused and compared for bulk water and the TEG water mixtures.The excess proton is found to reside predominantly on watermolecules but also resides on hydroxyl groups of TEG. Thelifetimes associated with structural diffusion time scales of theprotonated water were found to be on the order of 1 ps. Allthree functionals studied support the presolvation requirementfor structural diffusion. The highest level of theory shows a reduction in the free energy barrier for water water proton transfer inTEG water mixtures compared to bulk water. The effect of TEG shows no strong change in the kinetics for TEG watermixtures compared to bulk water for this same level of theory. The excess proton displays burst-rest behavior in the presence ofTEG, similar to that found in bulk water. We find that the TEG chain disrupts the hydrogen-bond network, causing the solvationshell around water to be populated by TEG chain groups instead of other waters, reducing the rigidity of the hydrogen-bondnetwork. Methylene is a dominant hydrogen bond donor for the protonated water in hydrogen-bond networks associated withproton transfer and structural diffusion. This is consistent with previous studies that have found the hydronium ion to beamphiphilic in nature and to have higher proton mobility at oil water interfaces.I. INTRODUCTIONPolyethylene glycol (PEG) (also known as poly(ethyleneoxide) (PEO) in its long chain form) is a water-soluble polymerthat has been widely studied for its interesting properties inaqueous environments. These properties include interactions ofits ethylene oxide (EO) chains in aqueous solutions,1 4 theability to adsorb on surfaces and interfaces,5 and protonconductivity in polymer electrolytes.6 18 The interaction ofwater with the PEG chain backbone has been widely studied todate, both by experiment2,3,19 and simulation4,20 26 due to theinteresting ability of PEG to form bridged hydrogen bonds(HBs) from EO monomer units to water molecules to otherEO monomer units within the same chain.4 These propertiesgive PEG a wide variety of applications such as biphasic (twophase) liquid liquid catalysis,27 electrochemical energy conversion as a proton exchange membrane additive,8 14,28 and indrug delivery via modification of therapeutic molecules (knownas PEGylation).29 31The conductivity of various ions with PEG in aqueoussolution has also been of research interest in past years. Inpharmaceutical applications, Capuano et el.32 studied thetransport properties of NaCl in PEG-water mixtures to 2016 American Chemical Societyunderstand the combined effect of PEG and salt interactionsin aqueous solution for extension to PEGylated proteins fordrug delivery. Consta and Chung33,34 used hybrid Monte Carloand classical molecular dynamics (MD) simulations todetermine charge-induced conformational changes and ionrelease mechanisms of PEG chains in aqueous chargednanodroplets for applications in atmospheric aerosols andelectrospray mass spectroscopy experiments. In fuel cellapplications, Honma and co-workers15 17 synthesized organic/inorganic nanocomposite hybrid polymer electrolyte membranes of SiO2/PEO via sol gel processes. Sundholm and coworkers35,36 synthesized and characterized PEO sulfonic acidsas a polymer electrolyte material. Ennari and co-workers37 40performed computational studies of acidic aqueous PEO andPEO sulfonic acid solutions of the similar molecular weight asthose of Sundholm and co-workers35,36 Classical MDsimulations were used to study the structure, solvation andproton diffusion to determine the atomic-level effect of theReceived: March 8, 2016Revised: May 24, 2016Published: May 24, 20165223DOI: 10.1021/acs.jpcb.6b02445J. Phys. Chem. B 2016, 120, 5223 5242
ArticleThe Journal of Physical Chemistry Bobserved that long-range proton transport occurs via periods of“burst” dynamics where the excess proton hops across multiplewater molecules for a larger distance on a short time scale,followed by periods of long-lived “rest” ( 5 ps) where it isstabilized by the local solvation structure on the picosecond(ps) time scale, still undergoing local proton hopping. Tofurther clarify the role of both the weak hydrogen bondformation to the hydronium ion and proton hopping overmultiple waters simultaneously to the molecular mechanismthat gives rise to long-range proton transport, Tse et al.62showed that concerted hops over multiple water moleculesoccur during both burst and rest periods and that the weakhydrogen bond formation to the hydronium was the strongercorrelation to the burst behavior. Thus, long-range protontransport occurs via single and multiple hops and is mostlycoupled to a favorable hydrogen bond network paired with aweak hydrogen bond formation to the hydronium to initiate theforward hopping, burst dynamics. From this refined state ofunderstanding of proton transport in bulk water, a naturalextension is to also study and refine proton transport inaqueous solutions and mixtures with other relevant molecularspecies that can have an effect on proton transfer mechanismsand transport.With advances in computational resources and algorithmsrelevant to electronic structure calculations, larger systems suchas study of proton transfer in aqueous polymer solutions can bemodeled via CPMD simulations. The work by Ennari and coworkers37,39,68,69 used classical MD simulations with anempirical force field that does not explicitly contain theelectronic degrees of freedom in its definition. This implies thatbond breaking and formation required for structural diffusion ofproton transport could not be simulated in such a way that isavailable in CPMD simulations. To overcome this limitation,the systems simulated included both a hydronium ion to modelthe vehicular diffusion and a particle with the mass and chargeof a proton that, interacting via electrostatics, could move freelyfrom surrounding water molecules to model the structuraldiffusion. This decomposition makes MD simulations of largersystems such as that required for long polymer chains viable butsuch simplification deviates from the real picture of structuraldiffusion. Mainly, it is not a single excess proton that hops fromone water molecule to another but a much larger, multidynamicand multitime scale molecular mechanism that is still a currentresearch challenge.For this study, we have simulated bulk water with one excessproton and a triethylene glycol (TEG) aqueous solution withone excess proton using CPMD simulations. This modelsystem provides further insight into the proton transport ofsystems with high PEG concentration and specifically the effectthat PEG has on proton transport compared to bulk water. Weanalyze the structures of the hydrogen bond networks that areformed, the different protonic defect species that occur andtheir associated lifetimes, solvation shell species and associatedlifetimes, mechanistic details of proton transfer, the free energybarriers for proton transfer between molecules, and thetransport kinetics and dynamics for proton transport.PEO chains and the sulfonic acid groups on the protonconductivity with varying water content. Yet, the theory of howproton transport occurs in aqueous solution began more thantwo centuries ago41 and the complete explanation of themechanism for proton transfer has only become fullyunderstood in the last two decades with the aid ofcomputational studies using a quantum treatment of thesystem.42 44Proton transport occurs via a combination of twomechanisms, the vehicular (center-of-mass translational motion) mechanism and the structural (proton hopping) diffusionmechanism, also known as the Grotthuss mechanism.45 48 Thelatter mechanism has been the elusive mechanism in understanding the anomalously high proton conductivity in aqueoussolutions. Common experimental techniques for studyingmicroscopic proton transport include nuclear magneticresonance45,46,48 and quasi-elastic neutron scattering,49 yetnanoscale-level resolution results beyond a qualitative level aredifficult to extract due to the atomic origin of protontransport.47 Yet, over the last two decades, computersimulations of liquids have given a quantitative understandingof the general hydrogen transfer reaction. Ab initio moleculardynamics (AIMD) simulations,50 specifically the “Car Parrinello” techniques for liquids,51,52 allow one to computethe forces on nuclei of atoms “on-the-fly” via electronicstructure calculations. This gives a technique that allows thestudy of systems free of parameter fitting, gives an atomic-leveldescription of the structure and dynamics as the system evolvesthrough time, and forces on nuclei that allow for the reactions(i.e., bond breaking and creation) to occur. This lastadvantageous quality of the Car Parrinello AIMD (CPMD)scheme is of central importance in proton transport in aqueoussolutions. Some of the pioneering CPMD simulations53,54determined that an excess proton in aqueous solution exists intwo previously postulated structures, the Zundel55 and Eigen56ion complexes, with the latter being the dominant species.Understanding the fundamental details of proton transportthrough aqueous medium continues to be a challengingresearch topic, both for experimental spectroscopic57 59 andcomputer simulation52,60 63 studies.Computational studies using AIMD have shown that protontransport is both a multidynamic and multitime scale processthat does not consist of a single proton hopping mechanism buta coupling of multiple molecular mechanisms.60,62,64,65 Thechemical bonding formation is coupled to the dynamicrearrangement of the hydrogen bonding network, which iscoupled to larger-scale motions that determine conditions forlong-range proton transport.62 In the work of Berkelbach etal.,60 the structural diffusion mechanism for proton transport isstudied by observing the change in solvation structure of thehydronium and the proton-receiving water,54,66 only for protonhops that then transfer to a different water molecule, thus notundergoing a reverse proton transfer reaction between theoxygen atoms. This stepwise, concerted mechanism for longrange proton transport is described as a weak hydrogen bondformation mechanism, where the hydronium begins to accept aweak, donating hydrogen bond from a water molecule while,simultaneously, the proton-receiving water loses an acceptinghydrogen bond, giving an exchange of solvation structurebetween the hydronium and the water. This exchange ofsolvation structure between the hydronium and the protonreceiving water is referred to as the “presolvation” concept.47,60,63,67 In a recent work of Hassanali et al.,65 theyII. SIMULATION DETAILSII.A. Simulation Size and Cell Setups. Density functionaltheory70 (DFT) based AIMD simulations, specifically CPMDsimulations, were performed in the microcanonical ensemble(NVE) using the Quantum ESPRESSO (opEn-Source Packagefor Research in Electronic Structure, Simulation, and5224DOI: 10.1021/acs.jpcb.6b02445J. Phys. Chem. B 2016, 120, 5223 5242
ArticleThe Journal of Physical Chemistry BOptimization) software suite71 for an excess proton in bulkwater and a TEG water mixture. For the bulk watersimulations, the system consisted of 32 water molecules andone excess proton in a cubic simulation cell of length 9.8652 Å,taken from previous studies.42 The TEG water mixtureconsisted of 27 water molecules, an excess proton, and oneTEG chain in a cubic simulation cell of length 10.0 Å. Bothsystems were subject to periodic boundary conditions for bothclassical and ab initio MD simulations.The system size was chosen to ensure strong interaction andobservable effect of the TEG chain with the excess proton. Fora larger system, the aqueous TEG solution becomes moredilute, decreasing the observation of the TEG effect on protontransport. Thus, our system size helps increase interactionsampling of possible proton transport effects induced by theTEG chain. The simulation cell sizes were determined afterclassical MD simulations were used for rough equilibration andgeneration of the initial nuclear coordinates in the system, to bediscussed further below. Using this system cell size, multipleAIMD simulations with different initial nuclear coordinates,were carried out to sample different phase space trajectories inorder to increase the probability of seeing TEG interactionswith the excess proton and also to increase the statisticalsampling of structure and transport of the excess proton.II.B. Classical MD for Generating Initial Conditions forAIMD. Classical MD simulations using the OptimizedPotentials for Liquid Simulations - All Atom (OPLS-AA)empirical force field72 74 for the TEG molecule and theTransferable Intermolecular Potential - Three Point (TIP3P)water model75 with flexible OH bonds76 for the water moleculeand the hydronium molecules were performed prior to theAIMD simulations to generate initial nuclear coordinates forinitialization of the AIMD systems. For the classical MDsimulations, the hydronium charges were taken from the workof Urata et al.77 The long-range Coulombic interactions werecomputed using the k-space style, long-range particle particleparticle-mesh solver78 with a relative error in the force toleranceof 10 4 kcal/mol. Geometric mixing rules were used forcomputing pair-style force field parameters of different atomtypes. The simulation cell was 1 charged due to thehydronium ion (not a charge neutral system). The pairwiseand Coulombic forces used a cutoff of 10.0 Å. To initialize thesystem, the TEG chain, water molecule and hydronium ionwere put on a simple cubic lattice arrangement of very lowdensity ( 0.3 g/cm3) initially to reduce aphysical initialconfigurations of high energy using the open-sourceMoltemplate molecule building and force field databasecode.79 The Large-Scale Atomistic/Molecular Massively ParallelSimulator (LAMMPS) open-source code was used to performall classical MD simulations.80,81 The reversible referencesystem propagator algorithm (rRESPA) by Tuckerman et al.82was used for integrating the equations of motion. The largesttime step was 1 fs (fs) for pairwise interactions and a smallesttime step was 0.2 fs for bonding interactions. The system wasinitialized by first performing energy minimization with anenergy convergence tolerance of 10 5 kcal/mol and a forceconvergence tolerance of 10 7 kcal/mol. After minimization,the system was run in the isobaric isothermal ensemble(NPT) for 250 ps using the Nosé Hoover thermostat andbarostat83,84 with the equations of motion formalism from thework of Shinoda et al.,85 the hydrostatic equations of Martynaet al.,86 the strain energy proposed by Parrinello and Rahman,87and the time integrator closely following the one fromTuckerman et al.88 The thermostat was set at 300 K with arelaxation timelength of 0.1 ps. The barostat was set at 1 atmwith a relaxation timelength of 1 ps with a drag applied to theNosé Hoover equations to dampen the large pressureoscillations that can occur in going from low density to highdensity. This protocol gave a converged density of 1.05 g/cm3for the system. This density appears to be within reasonableagreement with experiment, lying between the density of purewater (0.99 g/cm3) and the density of pure TEG (1.12 g/cm3)at STP. With a TEG mole fraction of 0.0357 and at 300 K,we are in relatively good agreement with the recent work ofBegum et al.,89 who report a density of 1.0412 g/cm3 for a molefraction of TEG equal to 0.05023 and at 303.15K for TEG water mixtures. Keeping the density fixed, the system thenunderwent thermal annealing via temperature ramping from300 to 900 K within 250 ps and then run in the canonicalensemble for a fixed temperature of 900 K for 500 ps and thenquenched back to 300 K for 500 ps. This allowed the system tocross possible high energy barriers that may not be overcome atlow temperatures and to explore more of the phase spacewithin the simulation time. The initial configurations for AIMDsimulations were taken from the final 500 ps at 300 K of theclassical MD simulations at one of the five equally spaced (100ps) configurations.II.C. AIMD. For the TEG water mixture, three differentDFT conditions were set up, each with either a differentexchange-correlation (XC) functional used in the quantummechanical treatment of the electronic degrees of freedom, adifferent pseudopotential used to describe the interactionbetween the valence electrons and the ions (nuclei and theircore electrons), or the inclusion of nonlocal electroncorrelation effects that are responsible for the van der Waals(vdW)/dispersion interactions, or vdW-inclusion/dispersionDFT. The first set of conditions, deemed the PW91 set, usedthe Perdew Wang 91 (PW91), generalized-gradient approximation (GGA) XC functional,90,91 Vanderbilt ultrasoftpseudopotentials for core electron-valence electron interactions,92,93 and no vdW/dispersion interactions. The second setof conditions, deemed the PBE set, used the Perdew, Burke,and Ernzerhof (PBE) GGA-type XC functional,94,95 Martins Troullier (MT) norm-conserving pseudopotentials,96 and novdW/dispersion interactions. The third set of conditions,deemed the PBE TS-vdW set, used were the same as the PBEset but included the Tkatchenko Scheffler density-dependentvdW/dispersion functional (TS-vdW).97,98 For the bulk watersimulations, four sets of DFT XC functionals were used: thePW91 XC functional with Vanderbilt ultrasoft pseudopotentials, the Beck, Lee, Yang, Parr (BLYP)99,100 XC functional withMT pseudopotentials, the PBE XC functional with MTpseudopotentials, and the PBE XC functional with theTkatchenko Scheffler density-dependent vdW/dispersionfunctional (TS-vdW).97,98 The choice of including the BLYPfunctional for the bulk water simulations was made to allow forcomparison with results in the literature.101,102For each set of DFT conditions, four to five independentAIMD simulations of the TEG water mixture were run,totaling 13 different simulations across all XC functionals. Sincethe system is positively charged due to the hydronium ion, anegatively charged jellium/background charge103 was applied tothe system to give a charge-neutral cell. The Kohn Shamorbitals were expanded in plane waves with a kinetic energycutoff 80 Ry for all three sets of DFT conditions. This cutoff isalso consistent with previous work for determining accurate5225DOI: 10.1021/acs.jpcb.6b02445J. Phys. Chem. B 2016, 120, 5223 5242
ArticleThe Journal of Physical Chemistry BFigure 1. Radial distribution functions of oxygen and hydrogen of water and hydronium ion in bulk water and TEG water mixtures at 300 K. Thedifferent figures are as follows: (a) O (water) O (water), (b) O (water) H (all hydrogens), (c) O (hydronium) O (water), and (d) O(hydronium) H (all hydrogens).cutoffs for pure water CPMD simulations.104 107 The chargedensity was expanded in plane waves with an energy cutoff of320 Ry for the PBE and vdW-PBE setups and 960 Ry cutoff forthe PW91 setup since ultrasoft pseudopotentials trade-off asmaller Kohn Sham orbital cutoff for a larger charge densitycutoff. The electronic subsystem was sampled at the Γ-point atthe first Brillouin zone. This is also consistent with previousliterature, specifically in the work of Galli et al.104,105 whereelectronic structure calculations of both a 64 water moleculecell and, by repeating the 64 molecule cell periodically twice ineach direction, a 512 water molecule cell, showed energy permolecule changes were less than 0.01 kcal/mol and the forcewas equivalent per atom within 10 5 a.u. between both systems.For the all the CPMD simulations, a time step of 0.075 fs anda fictitious electronic mass of 300 au was chosen due to testingfrom previous studies104 107 of CPMD for water simulations.The small time step and electronic mass are due to the lighthydrogen nucleus within the system. To ensure adiabaticitybetween the ions and the electrons, the fictitious kinetic energyof the orbitals was monitored and showed to be constant andremained close to zero. This test is a measure of how close thesystem stays to the ground-state Born Oppenheimer (BO)energy surface.50 If incorrect parameters are used, the electronswill exchange energy with the ions (destroying adiabaticity) andthis kinetic energy will increase, undermining the integrity ofthe CPMD simulation.For the bulk water AIMD simulations, the simulations werefirst run with 10 ps in the canonical (NVT) ensemble andfollowed by 50 ps of NVE for data production for each XCfunctional. The TEG water mixture AIMD simulations werefirst carried out in the NVT ensemble for 13 ps of equilibrationand then run in the NVE ensemble for a minimum of 10 ps fordata production for each of the simulations. The length of eachsimulation varied but it was ensured that 50 ps of dataproduction per XC functional set were obtained. The collectivetotal for each XC functional were 95 ps for the PBE set, 68ps for the PW91 set, and 52 ps for the PBE TS-vdW set. Thesystem was initialized by first relaxing the electronic wavefunction by relaxing the electronic degrees of freedom withfrozen, immobile ions. Then, the system wave function wasrelaxed by mobilizing the ions, allowing the ionic and electronicdynamical degrees of freedom to converge on the ground-stateenergy simultaneously. This was performed via geometryoptimization using dampened dynamics with a total energyconvergence tolerance and a fictitious electronic kinetic energyconvergence tolerance of 10 7 and 10 4 Hartree a.u.,respectively. For the equilibration stage of the TEG watermixtures in the NVT ensemble, Nosé Hoover chain thermostats were implemented,108 each with a chain length of four. Toachieve rapid equipartition of the thermal energy, one Nosé Hoover chain was applied per atom (“massive” Nosé Hooverthermostat). Four different thermostat frequencies were used tostimulate both the slow and fast intramolecular, oxygen hydrogen vibrational motions. The ionic temperature wasmonitored and was found to be well-averaged around 300 Kafter 3 ps. This was the followed by an additional 10 ps of NVTequilibration and then switched to the NVE ensemble for aminimum of 10 ps for data production.III. RESULTS AND DISCUSSIONIII.A. Comparative Analysis of Solvation Structure viaRadial Distribution Functions. Our analysis begins with adescription of the solvation structure found around H2O andH3O molecules for both the bulk water and TEG watermixture systems using radial distribution functions (RDF), g(r).We present a subset of all RDFs calculated from the simulations5226DOI: 10.1021/acs.jpcb.6b02445J. Phys. Chem. B 2016, 120, 5223 5242
ArticleThe Journal of Physical Chemistry BFigure 2. Radial distribution functions of oxygen for water and hydronium ion in TEG water mixtures with the oxygens of TEG chain at 300 K. Thedifferent figures are as follows: (a) O (water) O (ether of TEG), (b) O (water) O (hydroxyl of TEG), (c) O (hydronium) O (ether of TEG),and (d) O (hydronium) O (hydroxyl of TEG).signifying that H atoms are less likely to be locked intohydrogen bonds with oxygen of water molecules. Also, we seethat inclusion of vdW/dispersion interactions also reduce thepeak height, indicating a further softening of the waterstructure.107 It has become well-known that the PBE level oftheory gives an overstructured oxygen oxygen RDF for waterdue to the key limitations in GGA-DFT. These limitationsmanifest in theoretical descriptions of water as overstructuring,giving a representation of supercooled liquid water at ambientconditions. Recently, DiStasio et al. determined that includingvdW/dispersion interactions strengthen the interaction between a given water molecule and those found in its first andsecond solvation shells.107 The stronger interaction with secondsolvation shell water molecules allow them to move inward tointerstitial regions, weakening the HBs in the first solvationshell and allowing first solvation shell water molecules to moveoutward. Thus, we attribute this softening of structure with thedispersion-included PBE TS-vdW level of theory due toweakened HBs and stronger interaction with second solvationshell species, displacing first solvation shell water moleculescompared to the PBE level of theory. The most evident featureis the increase in the RDF value for the first minimum in thegO(H2O) O(H2O)(r) in Figure 1a. This is due to interstitial watermolecules between the first and second solvation shell, i.e. thewater molecules migration outward, from first to secondsolvation shells, and inward, from second to first solvationshells.For the solvation structure around H3O , we analyze theRDF of O in H3O with O in H2O, gO(H3O ) O(H2O)(r) (Figure1c), and with H atoms, gO(H3O ) H(r) (Figure 1d). IngO(H3O ) O(H2O)(r), we observe the first peak around 2.5 Å forboth systems and all levels of theory. The first minimum occursto illustrate key differences in local liquid structure between thebulk water system and the TEG water mixture, and also tohighlight affinity of H2O and H3O for specific groups along theTEG polymer chain specific to the TEG water mixture. InFigure 1, the RDFs for oxygen of H2O and H3O with differentspecies are shown.For understanding the solvation structure around watermolecules, we first analyze the correlation between the O ofH2O, annotated O (H2O), with other O (H2O) RDF, orgO(H2O) O(H2O)(r) in Figure 1a. The oxygen oxygen RDF inFigure 1a shows that the first peak around 2.75 Å is the samefor the bulk water and the TEG water mixture and all densityfunctionals. The first minimum occurs around 3.2 3.3 Å withthe location being found at greater distance from the oxygenwith an increasing level of theory for the XC functional. Thesecond peak occurs at around 4.2 4.3 Å, appears to vary lessfor the XC functional in the TEG water mixture.In Figure 1b, the RDF between the oxygen of water, O(H2O), and surrounding H, gO(H2O) H(r) is shown. The firstpeak occurs due to covalently bonded hydrogen atoms. We seethat the second and third peaks occur around 1.75 and 3.2 Å,respectively, for both bulk and TEG water mixture, regardlessof the XC functional. The peak height decreases with increasinglevel of theory and from bulk water to the TEG water mixture,indicating a weakening of the HB network.Based on the RDFs for oxygen atoms of water, the presenceof TEG shows an overall reduction in RDF peak height andsoftening of the local water structure. This occurs most likelydue to the TEG chain species being located in the first solvationshell, causing softening of the first solvation structure anddisrupting the HB network. This disruption is best shown inthe reduction of the second peak of the O of H2O and any H,5227DOI: 10.1021/acs.jpcb.6b02445J. Phys. Chem. B 2016, 120, 5223 5242
ArticleThe Journal of Physical Chemistry BTable 1. Protonic Defect Species Occurrences and Lifetimes Found in Aqueous H and Triethylene Glycol (TEG) SystemsafunctionalPW91proton rattlingprotonic defect%Σ%lifetime (fs)includedH3O ROH2 bOH RO cH3O ROH2 bH3O ROH2 bOH RO cH3O ROH2 bH3O ROH2 bOH H3O ROH2 .1399.970.03 9.9910099.8710099.9710010099.9710034.6 ( 55.5)34.6 ( 39.4)15.6 ( 13.4)7.6 ( 3.9)700 ( 1380)475d47.6 ( 94.0)42.7 ( 96.2)14.8 ( 12.77)7.6 ( 3.9)1290 ( 2430)110d45.9 ( 81.1)8.2 ( 1.5)1.6d1300 ( 2100)126.61dexcludedPBEincludedexcludedPBE TS-vdWincludedexcludedaProton rattling effects are included and excluded. bProtonated hydroxyl group of the TEG chain. cDeprotonated hydroxyl group of the TEG chain.Single lifetime measurement.dseem to have a drastic decrease in RDF size, going from thePBE to the PBE TS-vdW XC functionals, indicating theimportance of including dispersion for this system. This showsoverall there is a softening of the second solvation shellstructure around H2O via including dispersion.For Figure 2c and d, we show the RDF for O of H3O withthe ether O and hydroxyl O of TEG, respectively. In Figure 2c,the ether O of TEG is shown to be in the first solvation shell ofH3O for the PBE TS-vdW XC functional but not for thePW91 and PBE levels of theory. There appears to be a stronginteraction between the ether oxygen of TEG and H3O that isnot captured by modeling without dispersion effects. In Figure2d, the hydroxyl O of TEG is in the first solvation shell of H3O for the PW91 and PBE levels of theory but not for the PBE TS-vdW XC functional. Interestingly, the PBE TS-vdW XCfunctional gives an almost feature-less RDF between a H3O and the hydroxyl group of TEG, whereas, the PW91 and PBERDFs have apparent peaks at 2.5 and 4.5 Å. The strongestinteraction of H3O with the hydroxyl group occurs with thePW91 level of theory. Overall, for H3O , the PW91 and PBEXC functionals give a H3O solvated mostly by H2O but also byhydroxyl O of TEG in the first solvation shell and PBE TSvdW gives H3O solvated by H2O and ether groups of TEG.III.B. Protonic Defect Species Occurrence and LifetimeAnalysis. One interesting effect that TEG has on protontransfer in aqueous solution is the different protonated anddeprotonated species, or protonic defect species, that can occurand their associated lifetimes. We analyze the trajectory forthese species by, first, assigning every hydrogen atom to thegeometrical
of the general hydrogen transfer reaction. Ab initio molecular dynamics (AIMD) simulations,50 specifically the "Car Parrinello" techniques for liquids,51,52 allow one to compute the forces on nuclei of atoms "on-the-fly" via electronic structure calculations. This gives a technique that allows the
