Article pubs.acs.org/JPCC
Probing the Structures of Hydrated Nafion in Different Morphologies Using Temperature-Accelerated Molecular Dynamics Simulations Jeremy Lucid,† Simone Meloni,*,†,§ Donal MacKernan,†,∥ Eckhard Spohr,‡ and Giovanni Ciccotti†,⊥ †
School of Physics, University College Dublin, Dublin, Ireland Lehrstuhl fuer Theoretische Chemie, Fakultaet fuer Chemie, Universitaet Duisburg-Essen, Germany
‡
ABSTRACT: We perform combined temperature-accelerated and standard molecular dynamics (MD) simulations to elucidate the atomistic structure of hydrated Nafion (hydration level λ = 6.5) in the slab and cylinder morphologies. Our samples are initially made of elongated Nafion strands with a relatively small fraction of gauche defects. Our simulations show that even very long (>50 ns) “brute force” MD simulations are insufficient to reach equilibrated structures. In fact, ∼30−40 ns long temperature-accelerated molecular dynamics (TAMD) simulations started from the same initial conditions explore more stable (lower potential energy) stationary structures. The effect of TAMD is to allow a rearrangement of the backbone consisting of an increase in gauche defects, which cannot be obtained by “brute force” MD because the trans−gauche transition is a rare event at room temperature. Associated with the backbone rearrangement, we observe a change in the structure of the water layers/tubes as measured by the size and number of bulk (four-fold coordinated water molecules) and surface-like water clusters. At equilibrium, the mean size of bulk-like water clusters is small, typically between 10 and 20 molecules, depending on the morphology. Larger clusters are also present in our samples, the largest being made of ∼350 molecules, but even the latter is too small for percolation. This suggests that the proton transport through each morphology might be a two-step process: Grotthuss-like within bulk-like water clusters and of a different type (e.g., diffusive or even transport across fluctuatively opening necks) between clusters.
1. INTRODUCTION
membrane via the same proton transport mechanisms active in water and aqueous acidic solutions. Thus, a key requirement for fuel cell operation is the membrane’s ability to provide adequate morphological features for efficient proton transport. These consist of an interconnected and percolating network of aqueous clusters (i.e., a network of clusters connecting two ends of a membrane subdomain), channels, and lamellar regions that internally exhibit water environments locally not too different from bulk aqueous systems. Ideally, such a network has an ordered structure that allows for fast proton transport through the membrane. In addition, the membrane must have the necessary structural integrity to permanently separate the anode and cathode region, which rules out many choices of membrane material. Polymer electrolytes that can provide these features and, at the same time, are sufficiently robust, both chemically, thermally, and mechanically, usually consist of sulfonated perfluorinated compounds, although other materials have been employed as well.2−6 The most widely employed chemical compound in the class of these perfluorinated sulfonic acids is Nafion by DuPont. Proton conduction is based on dissociation of protons from the highly acidic sulfonic acid (−SO3H) groups when brought in
Hydrogen/oxygen (H2/O2) or polymer electrolyte fuel cells (PEFCs) as well as direct methanol fuel cells (DMFCs) are electrochemical energy conversion devices that combine fuel oxidation at the anode and oxygen reduction at the cathode to produce electricity from the chemical energy stored in the fuel. While energetically downhill (and therefore spontaneous), the elementary reaction steps are not fast enough to achieve adequate power output without the use of catalysts. Consequently, the elementary electrocatalytic processes at the anode and cathode are presently the focus of extensive research using a multitude of thermodynamic, spectroscopic, electrochemical and scanning probe experiments, and computational techniques based on quantum mechanical density functional approaches.1 The long-term goal of this research is to replace the use of precious metal catalysts with cheaper materials. A second requirement for this electrochemical process is the spatial separation of oxidation and reduction reactions, achieved through the separation of the anode and cathode by a membrane. Ideally, such a membrane is permeable only for one ionic species. In PEFCs and DMFCs this membrane is an electrolyte consisting of a polymer with chemically attached anionic groups and dissociable protons. While the anions (usually sulfonate groups, SO3−) are stationary during fuel cell operation, the protons, when brought in contact with water, can rapidly be transported through the © 2012 American Chemical Society
Received: September 11, 2012 Revised: December 15, 2012 Published: December 18, 2012 774
dx.doi.org/10.1021/jp309038n | J. Phys. Chem. C 2013, 117, 774−782
The Journal of Physical Chemistry C
Article
contact with water. The remaining sulfonate groups SO3− are chemically bound at the end of short side chains attached to a polymeric backbone (see Figure 1). The current understanding
Before one is able to address the morphological changes resulting from the coupling between flow and local pore structure, one first needs to understand the morphology of Nafion membranes in equilibrium. In the past, many attempts have been made in this direction, both within the framework of atomistic MD simulations and, recently, dissipative particle dynamics (DPD) simulations33,34 with ever increasing system sizes, which, however, still do not exceed the 100 nm scale. Given that a Nafion 117 (see Section 2.1) membrane has a thickness of 175 μm, it is clear that only small models of a real fuel cell membrane can be studied. An additional complication is the disordered nature of the polymer electrolyte, and the extremely long relaxation times of the polymer network as compared with those of the local proton transfer processes. To date, MD studies of proton transport and the morphology of aqueous domains have investigated only local (meta)stable domains that did not exhibit any substantial changes in simulations of some tens of nanoseconds in duration.35−58 Certainly, without further examination of such computer-generated configurations one cannot assume that the generated morphologies are even marginally stable and relevant in the statistical-mechanical sense. Before studying morphologies in steady-state flux or time-dependent situations it is important to assess the morphological stability of equilibrium structures. The purpose of the present exploratory work is to investigate this issue of stability and structure of local membrane morphologies by using temperature-accelerated molecular dynamics59 (TAMD), a simulation technique allowing an efficient exploration of the configuration space of a system in the presence of metastabilities. In principle, TAMD allows the reconstruction of the (Landau) free energy of a set of suitable collective coordinates. One can then identify the minima of this free energy and characterize the atomistic structure of configurations associated to them. However, despite recent improvements,60 reconstructing with TAMD the free energy of very many collective variables (CVs) (thousands, in the present investigation) is computationally inefficient. Thus, we use TAMD only to explore the configuration space and then use “standard” MD, started from a low potential energy configuration identified in the first step, to characterize the properties of highly probable states. The remainder of this article is organized as follows. In Section 2, we detail how the initial structures of the atomistic samples were generated (Section 2.1) and describe the computational methods used in the present investigation (Sections 2.2, 2.3, and 2.4). In Section 2.5, we report the details of our simulations and in Section 3 present and discuss our results on the structure of the slab and cylindrical morphologies (Section 3.1) and those of the associated water layers/tubes (Section 3.2). Finally, in Section 4, we draw some conclusions and give an outlook for future work.
Figure 1. Structure of the monomeric unit forming Nafion. In the present work, l = 6, m = 1, and n = 1. The polymer chains present in our samples contain 20 monomers.
of proton transport in Nafion is based on extensive evidence from small-angle neutron and X-ray scattering experiments and is supported by atomistic and mesoscopic computer simulations (see, e.g., refs 2 and 7−16 for reviews). Conducting protons are transported through aqueous clusters (ponds), channels, or (locally) slab-like regions, which are only established in the presence of a sufficient amount of water. When the water content in these ponds becomes high enough, these structures become interconnected and proton transport is sufficiently fast to generate high power output from a fuel cell. A large number of mostly molecular dynamics (MD) simulation studies of proton transport in a variety of geometrical models of local morphology have been performed in the past, using classical - both nonreactive and reactive - force field approaches or, more recently, ab initio MD simulations of small channel fragments.17−23 In the case of proton exchange membranes, the design goal is to allow for proton transport without additional water transport, which would only be possible if classical ion transport could be totally suppressed. However, proton transport is intrinsically connected to water mobility via the coupling of classical ion transport (in which the proton is transported as a hydronium ion, H3O+, Zundel ion, H5O2+, or Eigen species, (H3O+)·3H2O) and of the so-called Grotthuss hopping mechanism (where protons are shuttled from one water molecule to the next one, provided that the local geometry is adequately preconfigured through thermal fluctuations; see, e.g., refs 24 and 25 for reviews). In all of these cases, some water is inevitably moved through the membrane from the anode to the cathode side of a fuel cell together with the protons. This phenomenon of so-called electro-osmotic drag26,27 is a central aspect of fuel cell operation, subsumed under the term “water management” (see, e.g., refs 28−31 for neutron and X-ray radiography experiments on mesoscopic and macroscopic water distribution in membranes). In essence, the effect leads to an overall reduction of water concentration near the anode. Because local proton mobility increases with increasing water content, this results in a suboptimal proton transport through the membrane. At the cathode, in contrast, there is an excess of water due to the flow and to the cathodic reaction. The overall flow of water needs to be re-equilibrated either by back-diffusion through the membrane or through appropriate engineering measures external to the fuel cell. This macroscopic transport of water through the membrane can even lead to morphological changes during fuel cell operation (see, e.g., ref 32 for a recent study of possible fieldinduced morphological effects). Thus, the transport processes may lead to a steady-state morphology (during continuous fuel cell operation) different from the stationary one of the ideal membrane or to time-dependent processes between dynamically evolving water regions within the membrane during load changes (e.g., during startup or shutdown).
2. THEORETICAL METHODS 2.1. Atomistic Models of the Nafion Samples. Two morphological models of Nafion were selected, the slab model61 and the parallel cylinder model,62 both of which are considered idealized models of local structures present in real Nafion membranes.2 The samples used in this study are composed of Nafion, water, and hydronium molecules. Nafion is composed of four subunits (see Figure 1), where −[CF2− CF2]l− and −[CF−CF2]m− units compose the backbone of the polymer and the [O−CF2−CF-CF3]n and [O−CF2−CF2− SO3] units compose the side chains. The equivalent weight, that is, the ratio between the molecular weight and the number of sulfonate groups, determines the properties of Nafion. In this 775
dx.doi.org/10.1021/jp309038n | J. Phys. Chem. C 2013, 117, 774−782
The Journal of Physical Chemistry C
Article
replicated nine times in the z direction to form the full doublelayer slab system. Simulations at 1000 K for 1 ns were run to remove replication artifacts. Figure 2 shows a snapshot of the slab system, viewed along the y direction, after equilibration at P = 1.0 atm and T = 300 K for 1.5 ns. Our parallel cylinder model consists of a 2 × 2 grid of aligned polymer tubes filled with water, with the side chains pointing inside the water channels. Each tube is made of nine strands forming the lateral surface of the cylinder, giving a total of 36 oligomer strands for the entire system identical to that in the slab system. The cylinder system is filled with water to the same hydration level as that of the slab system, resulting in a water channel diameter of ∼1.4 nm with box dimensions in the x, y, and z directions of 5.52, 25.25, and 5.52 nm, respectively. The cylinders forming the system are constructed manually by placing elongated Nafion strands on the surface of a cylindrical scaffold, with the hydrophilic side chains pointing inward. No connecting bridge structures between adjacent cylinders were created in the initial sample. A top view along the y direction of the parallel cylinders system after equilibration at P = 1.0 atm and T = 300 K for 1 ns is shown in Figure 2. 2.2. Methods for Exploring the Configurational Space. In the present work, we combine TAMD64 with NVT MD59 to efficiently obtain statistical information over the relevant regions of the configuration space associated with the slab and cylinder models. In fact, MD is inefficient at sampling the configuration space when the system has metastabilities; that is, the free-energy landscape presents several minima separated by large free-energy barriers. This problem can be solved by employing TAMD, a method that allows exploration of the space of a set of CVs. However, especially when one uses a large number of CVs, obtaining accurate statistical information on the relevant part of this space might be inefficient (see below). This local exploration can be performed more efficiently by MD, which is why we combine the two techniques. In the following, we will briefly summarize the TAMD method. The TAMD method involves the introduction of an extended dynamical system consisting of the atoms together with a set of additional dynamical variables, denoted by z in the following. The z are the possible values of a set of CVs {θj(x)}j=1,n, which are suitable to characterize the system and the rare events associated with transitions between metastable states. This extended system evolves according to the following set of coupled equations
work, we focus on Nafion of 1100g/mol SO3, corresponding to commercial Nafion 117, which is the most commonly used membrane material in fuel cells. Nafion is a stochastic polymer; that is, the values of l and m in Figure 1 can change from molecule to molecule and, within the molecule, at different positions along the backbone, whereas n takes value one. In this work, we consider the simpler case of a regular polymeric structure in which two side chains are separated by six [CF2− CF2] units, that is, l = 6, and each Nafion polymer has 20 monomeric units. Such a molecule has an equivalent weight = 1046 g/mol SO3 close to that of Nafion 117. In this study, both model systems are constructed from elongated Nafion oligomers with a typical end-to-end distance of ∼25 nm. The choice of elongated as opposed to coiled polymer strands is an arbitrary one in that morphological scaffolds can, in principle, be constructed from either or both strand types. However, because it is unclear how system properties vary between coiled and elongated polymer scaffolds, we chose to confine our investigation to scaffolds constructed from comparable elongated strands to reduce the complexity of the problem. In this study, the slab and cylinder model scaffolds are made of an identical number of polymer strands, water content level, and comparable aspect ratio. This, together with the use of an analogous structure of the strands for both scaffolds, allows for a fair comparison of the characteristics of the two systems. In our slab sample, 36 Nafion strands are organized into two slabs initially separated by water layers of ∼1.1 nm in thickness, subject to periodic boundary conditions. The sample contains 3960 water molecules and 720 hydronium ions (H3O+), corresponding to a hydration level of λ = 6.5, where λ has its conventional definition of the number of water plus hydronium molecules per sulfonic acid group. Typical operational values for λ range from 2 to 16, where experimental studies63 have shown that the higher figure corresponds to a fully hydrated system, so that the water content of our sample is in the low water loading regime. Overall, the system is composed of 57 400 atoms with edge lengths in the x, y, and z directions of 5.18, 25.25, and 5.83 nm, respectively. The sample is created by first forming a slab “slice”, with dimensions 5.18, 25.25, and 0.65 nm in the x, y, and z directions, respectively. This consists of four equivalent parallel strands in pairs, aligned along the y axis with each pair being separated from the other pair of strands by a water layer (due to the periodic boundary conditions, our system is composed of two water layers, see Figure 2). This thin double-layer slab is then
⎧ ∂U (x , z) + thermo(β) ⎪ mixï = − k ∂xi ⎪ ⎨ ∂Uk(x , z) ⎪ + thermo(β ̅ ) ⎪ μj zj̈ = − ∂zj ⎩
(1)
where mi and xi denote the mass and position of the ith atom. μj is the “inertia” of the auxiliary variable zj, a tunable parameter, the effect of which will be described below. Uk(x,z) is the potential governing the system, consisting of the sum of two terms: V(x), the physical potential, and (1/2)k(θj(x) − zj)2, the term coupling x and z. Finally, thermo(β) and thermo(β̅) are the thermostats controlling the inverse temperatures of the atoms (β = 1/kBT) and of the additional degrees of freedom (β̅ = 1/kBT̅ ), respectively. Under the condition of adiabatic separation, that is, if z is much slower than x, which can be achieved by properly tuning μ and the thermostat parameters, and for k sufficiently large, the above dynamics samples the probability Pβ̅θ(z) = exp[−β̅F(z)], where
Figure 2. Two hydrated Nafion morphologies: (A) slab model and (B) cylinder model. The carbon atoms of Nafion are shown in green, sulfur in red, and oxygen atoms of water in blue. The densities of the slab and cylinder systems are 1.82 and 1.80 g/cm3, respectively. These densities have to be compared with the experimental value of 1.86 g/cm3 for samples at λ = 7.63 776
dx.doi.org/10.1021/jp309038n | J. Phys. Chem. C 2013, 117, 774−782
The Journal of Physical Chemistry C
Article
F(z), given below, is the Landau free energy at the physical temperature T n
F(z) = β −1 ln Z −1
∫ e−βV (x) ∏ δ(θj(xi) − zj) dx j=1
(2)
with Z = ∫ e−βV(x) dx. Even though the z are distributed according to the free energy F(z), they move on the free energy surface at the temperature T̅ . Thus, by setting T̅ ≫ T, the exploration of F(z) is more efficient because the system can more easily overcome the free energy barriers present in the system. One could, in principle, use TAMD to compute F(z) via the Boltzmann inversion of the histogram of the auxiliary variables z (F(z) = β̅−1 ln(Pβ̅θ(z))), but this presents two problems. First, Pβ̅θ(z) is more uniform than Pβθ(z), which implies that a much more accurate estimate is needed to obtain reliable values of F(z). Second, this approach scales exponentially with the number of CVs. Therefore, this procedure cannot be followed in a case, like the present one, in which a large number of CVs are used. Thus we use TAMD just to explore more extensively the conformational space. Then, the configuration space associated with a given region of the CV space is sampled by MD simulations starting from configurations explored by TAMD. 2.3. Collective Variables. Previous long simulations (∼30 ns) performed on analogous systems57 have shown that the rearrangement of the polymer side chains plays a central role in the relaxation of the proposed initial morphology of Nafion samples. For example, side chains were observed participating in the formation of water channels connecting adjacent cylinders within cylinder models and neighboring water layers within slab models. The rearrangements of Nafion side chains also lead to the roughening of the water clusters surface, quantified by a doubling of the solvent-accessible surface area.57 These observations suggest that observables monitoring the orientation of the side chains can be suitable CVs to explore the conformations of the system within the slab and cylinder morphologies. Thus, to each polymer strand we have associated 20 3D vectorial CVs corresponding to the center of mass of the side chains, resulting in a total of 720 3D CVs. 2.4. TAMD Parameters Setup. To perform TAMD calculations, we need to define the value of the inertia for the auxiliary variables, μ, the coupling constant, k, between the auxiliary variables and the CVs associated with the atomic system, and the parameters controlling the thermostats. As thermostats, we use Nosé−Hoover chains (with a chain length of three), and thus we need to set the inertia Qcv and Q of the chains coupled to the CVs and atoms, respectively. The value of k must be such that the estimated mean force, ▽FNk (z) = 1/NΣi=1,Nk·(θ(xi)−z), where xi is the configuration of the atoms at the ith MD step and N is the total number of steps, is well converged. A strategy for establishing a suitable value for k consists in computing the estimate of the mean force at selected z for different values of k. This can be obtained by performing MD simulations according to the first equation of motion of eq 1 at fixed z. In Figure 3, we report a few components of the running average of ▽FNk (z) (i.e., ▽FNk (z) vs N) at k = 50, 100, 150, 300, and 500 kcal/mol·Å2 for z = θ(x) corresponding to the initial configuration of our simulations. Notice, that the minimal k value for which the estimate of the mean force is well-converged is k = 150 kcal/mol·Å2 and so was the value set for the simulations presented in this manuscript. In fact, the lowest k value for which the mean force is well-converged is the optimal choice for this parameter because it allows for the longest possible time step as the
Figure 3. Mean force, ▽FNk (z), shown for three randomly chosen CVs, obtained for increasing values of the spring constant k.
biasing force is less stiff than with higher k values, which would require smaller timesteps to avoid numerical instabilities. Let us now move to setting the inertia of the auxiliary variables. μ must be chosen such that the CVs are much slower than the atoms. In practice, the atoms must be at the equilibrium for each given value of z. In this work, we exploit the fact that our CVs are the center of mass of the side chains, and thus their average speed can be obtained from the equipartition theorem: |ẋ| = ()1/2 = (2kBT/M)1/2, where ẋ is the velocity of the center of mass of a given side chain and M is the corresponding mass. TAMD requires that |ż| ≪ |ẋ|, and so we impose that |ż| = |ẋ|/30. The 1/30 ratio between the speed of the center of mass of the side chains and of the z variables ensures a good time scale separation between the two systems but does not make the sampling of the z space too slow as to be computationally unaffordable. We can now set the inertia of the z such that, given the temperature T̅ , their average speed is the |ż| = |ẋ|/30: μ = 2kBT̅ /|ż|2 = 302MT̅ /T Concerning the inertia of the Nosé−Hoover chains, we set Q to 9.7 × 108 kcal·fs2/mol, which is a typical value for systems of the kind considered here. With this choice, the mean force converges, typically, within 40 ps (see Figure 3). Thus, Q z is set such that the characteristic time of the CV thermostat (see Appendix B of ref 65) is 40 ps, corresponding to an inertia Q z = 1.7 × 1012 kcal·fs2/mol. 2.5. Simulation Details. All simulations were carried out using the fully atomistic model of Nafion 117, as described in Section 2. For both model systems, the following MD and TAMD simulations were performed using the LAMMPS (largescale atomic/molecular massively parallel simulator)66 code. The force-field parameters were chosen to be the same as those in ref 41, with the particle mesh (PPPM)67 method used for the electrostatic interactions. The equations of motion were integrated using the RESPA64 algorithm with an outermost time step of 2.0 fs and an innermost time step of 0.25 fs. Finally, atom indexes were reordered on-the-fly, assigning contiguous values to atoms belonging to the same linked cell.68 The entries of position, velocity, force, and so on vectors are reordered accordingly. After reordering, data of close in space atoms are also close in computer memory. This improves the efficiency of the force calculation and of the parallelization, resulting in a speed-up factor of 2 to 3. We began our investigation with pure NVT MD runs starting from the configurations described in Section 2. After a suitable equilibration time, we gathered equilibrium statistics for both. The information obtained from these runs is a benchmark to 777
dx.doi.org/10.1021/jp309038n | J. Phys. Chem. C 2013, 117, 774−782
The Journal of Physical Chemistry C
Article
Figure 4. Potential energy timeline of the (A) slab and (B) cylinder model along the various branches of our simulations.
which we can compare the TAMD results. Many short TAMD simulations were performed to find a suitable auxiliary temperature at which to study the model systems. Extensive TAMD simulations, of which we describe the results below, were then ran at the chosen β̅. We performed one extensive TAMD simulation per morphology, beginning from the same initial configuration of the MD simulations. Finally, we ran an additional MD simulation starting from a configuration extracted along the TAMD trajectory. This was done for the purpose of obtaining equilibrium statistics in the lowest potential energy configuration visited by TAMD. The three classes of simulations are listed and described in detail in the following. • MD from initial configuration (MD-IC): MD simulations in the NVT ensemble were performed at the average volume of the NPT simulations described in Section 2.1 for both slab and cylinder systems. For these simulations, the temperature of the atoms was maintained at 300 K. - Slab system (S-MD-IC): Total duration of 51 ns. Equilibrium statistics collection started at ∼42 ns, in correspondence with the achievement of stationarity of the potential energy. - Cylinder system (C-MD-IC): Total duration of 65 ns, with equilibrium statistics gathered after ∼50 ns.
- Cylinder system (C-MD-TC): MD simulations began after 32 ns of TD-IC and were performed for 12 ns.
3. RESULTS AND DISCUSSION 3.1. Potential Energy and Polymer Backbone Configurations. We are interested in investigating if TAMD simulations explore significantly different states from those readily visited along MD trajectories. To this end, we use the potential energy of the physical system for establishing whether it has reached the stationarity and whether TAMD explores more stable (lower potential energy) configurations. Figure 4 shows the timeline of the potential energy for the slab and cylinder systems during the simulations described in Section 2.5. The potential energy values and those of other observables, presented later, computed along MD and TAMD trajectories are consistently shown in blue and red, respectively. This Figure clearly shows that the MD-TC configurations are at a lower potential energy than the MD-IC configurations. For the slab system, S-MD-IC leads to a 7.6% decrease in potential energy from an initial value of −0.7 kcal/mol·atom to a final local time average value of −0.753 kcal/mol·atom, whereas S-MD-TC gives a final local time average value of −0.773 kcal/mol·atom, a 10.4% decrease in the potential energy. Similarly, for the cylinder system, which began at −0.65 kcal/ mol·atom, C-MD-IC gives an 8.9% decrease to a final local time average value of −0.708 kcal/mol·atom, whereas C-MD-TC results in a 17.1% decrease to a final average value of −0.761 kcal/mol·atom. In conclusion, the TAMD simulations performed using these CVs explore lower potential energy configurations than those obtained by standard brute-force MD simulations, while remaining within the desired morphology as determined by visual inspection. The strong oscillatory behavior shown by the potential energy curve along the C-TD-IC simulation is associated with an oscillatory bending motion of the four cylinders. This collective bending was also observed for the first 10 ns of S-TD-IC. These bending modes are consistent with previous observations,57 which show that idealized rigid and straight parallel cylindrical scaffolds evolve into slightly more bent tubular structures. To investigate the primary source of this decrease in potential energy along the TAMD trajectories, we focused our attention on two observables. First, we measured the mean distance between the sulfur and the backbone oxygen atoms of each side chain so as to test whether the drop in potential energy was associated with a more elongated or bent structure of the side chains. No correspondence was found between the
• TAMD from initial configuration (TD-IC): TAMD simulations began from the same initial condition as that of the MD-IC simulations. The temperature of the atoms was again kept at 300 K while the artificial temperature of the auxiliary system was set to 1000 K. The auxiliary temperature was chosen low enough so as to confine the exploration within the chosen morphology. Higher T̅ (3000 and 5000 K) was found to bring the systems to configurations deemed inconsistent with slablike or cylindrical-like morphologies. - Slab system (S-TD-IC): Total duration of 27 ns. - Cylinder system (C-TD-IC): Total duration of 38 ns. • MD from a TAMD configuration (MD-TC): MD simulations were performed beginning from one of the lowest potential energy configurations obtained along the TD-IC simulations. This allowed the physical system to relax and equilibrium statistics to be gathered. - Slab system (S-MD-TC): MD simulations began after 19 ns of TD-IC and were performed for 12 ns. 778
dx.doi.org/10.1021/jp309038n | J. Phys. Chem. C 2013, 117, 774−782
The Journal of Physical Chemistry C
Article
Figure 5. Instantaneous value of the mean 1−4 carbon distance for the (A) slab and (B) cylinder model.
great interest because this might suggest differences between proton transport mechanisms in the two systems. We analyzed water domains according to the amount of hydrogen-coordination related bulk-like (BL) or surface-like (SL) water molecules, together with the size and number of clusters formed by them. We define a BL water molecule as one forming four hydrogen coordinations with other water molecules, like in bulk water. This is different from the usual geometrical definition of a “Bulk” molecule. A molecule forms a hydrogen coordination if the oxygen−oxygen distance with a neighboring molecule is
