A molecular dynamics investigation ofwater and ion transport throughmodel carbon nanotubesbyLin LiuB.Sc., Nankai University, 2011A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFDOCTOR OF PHILOSOPHYinThe Faculty of Graduate and Postdoctoral Studies(Chemistry)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)March 2017c© Lin Liu 2017AbstractIn this dissertation, we investigate water and ion transport through carbon nanotubes usingmolecular dynamics simulations. Specifically, we examine how different water models influ-ence the simulated conduction rates. We consider three common water models, which areTIP4P/2005, SPC/E, and TIP3P, and observe that water flow rates through the same nan-otube are strikingly different amongst the different water models. Also, the water flow ratedependence on temperature fits an Arrhenius-type equation over a temperature range from260 to 320 K. We provide evidence that there are two factors which determine the conductionrate: the bulk fluid mobility, and the molecular structure of confined water. For narrow nan-otubes, for example, a (6,6) nanotube, where water only forms a single-file configuration, thefirst factor can largely account for the flow rate differences. In this case, we show that the con-duction rate correlates with the diffusion coefficient of bulk water. Our simulation results arewell described by continuum hydrodynamics as well. The factor of bulk fluid mobility is stillimportant in the water conduction through intermediate-size nanotubes, such as a (9,9) nan-otube. Also, the formation of complex configurations within such nanotubes can impede thetransport rate by influencing the mode of water conduction. The ordered structure occurringwithin nanotubes can also explain the differences between simulation results and continuumhydrodynamics predictions. Hence, both factors decide the water conduction rates throughintermediate-size nanotubes. Moreover, we demonstrate that the ion flow rate depends on theviscosity of the bulk solution, as well as the water structure within the nanotubes, togetherwith the ion size. In particular, at lower temperatures complex water configurations act toimpede ion transport while still allowing water to flow at a significant rate. In general, ourefforts on this issue are of importance for future simulation studies investigating water andiiAbstraction conduction through nanoscopic channels. This dissertation might also prove useful indesigning more efficient nanoscopic conduits for future experimental studies.iiiPrefaceThe research presented in this dissertation has appeared as co-authored, peer-reviewed pub-lications by L. Liu and G. N. Patey. Parts of this dissertation have been published as journalarticles:• L. Liu and G. N. Patey, “Simulations of water transport through carbon nanotubes: Howdifferent water models influence the conduction rate”, J. Chem. Phys., 141, 18C518(2014).• L. Liu and G. N. Patey, “Simulated conduction rates of water through a (6,6) carbonnanotube strongly depend on bulk properties of the model employed”, J. Chem. Phys.,144, 184502 (2016).• L. Liu and G. N. Patey, “A molecular dynamics investigation of the influence of waterstructure on ion conduction through a carbon nanotube”, J. Chem. Phys., 146, 074502(2017).The computer simulations and data analysis were conducted by L. Liu with guidance andsuggestions from G. N. Patey. The manuscripts were written by L. Liu with revisions andpolishing by G. N. Patey.ivTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiiList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xList of Symbols and Abbreviations . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xv1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Molecule Transport through CNTs . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Computational Studies of Molecule Transport through CNTs . . . . . . . . . . 31.3 Water Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61.4 Outline of Dissertation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 Models and Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102.1 Modelling of Water . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102.2 Modelling of CNTs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112.3 Interaction Potentials . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132.4 Simulation Details . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14vTable of Contents2.5 Application of Pressure Difference . . . . . . . . . . . . . . . . . . . . . . . . . 172.6 Data Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 182.6.1 Flow Rates and Entry Rates . . . . . . . . . . . . . . . . . . . . . . . . 192.6.2 Characterization of Confined Configurations . . . . . . . . . . . . . . . 192.6.3 Diffusion Coefficient . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 202.6.4 Potential of Mean Force . . . . . . . . . . . . . . . . . . . . . . . . . . . 213 Simulated Water Conduction Rates through Intermediate-size CNTs . . . 223.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 223.2 Water Flow Rates through Intermediate-size CNTs . . . . . . . . . . . . . . . 223.3 Water Structure inside Intermediate-size CNTs . . . . . . . . . . . . . . . . . . 233.4 Water Transport Mechanisms . . . . . . . . . . . . . . . . . . . . . . . . . . . . 313.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 354 How Water Conduction Rates through CNTs Is Related to Bulk Propertiesfor Different Water Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 364.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 364.2 Water Flow Rates through a (6,6) CNT . . . . . . . . . . . . . . . . . . . . . . 374.3 Explanation of the Model-dependent Water Transport Rates . . . . . . . . . . 394.4 Water Hydrodynamics and Bulk Fluid Properties . . . . . . . . . . . . . . . . 434.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 545 How Different Water Models Affect Simulated Ion Transport Rates througha (9,9) CNT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 565.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 565.2 Water and Ion Flow Rates through a (9,9) CNT . . . . . . . . . . . . . . . . . 575.3 Why Water and Ion Conduction Is Dependent on the Water Model Employed 615.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72viTable of Contents6 Conclusions and Perspective . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 746.1 Simulation Results of Water and Ion Transport through CNTs . . . . . . . . . 746.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79AppendicesA Molecular Dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90A.1 Ewald Summation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90A.2 System Evolution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92A.3 Thermostat . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93A.4 Potential of Mean Force . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94B Discussion of Some Simulation Details . . . . . . . . . . . . . . . . . . . . . . 97B.1 Cell Dimension . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97B.2 Thermostat . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97B.3 Implementation of Pressure Difference . . . . . . . . . . . . . . . . . . . . . . . 98viiList of Tables2.1 List of geometric parameters of water models . . . . . . . . . . . . . . . . . . . 102.2 List of force field parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153.1 Water flow rates through (8,8) and (9,9) nanotubes . . . . . . . . . . . . . . . . 233.2 Percentage of ring-bound/ring-free water in (8,8) and (9,9) nanotubes . . . . . 273.3 Numbers of hydrogen bonds per water in (8,8) and (9,9) nanotubes . . . . . . . 314.1 Water flow rates through/entry rates into nanotubes and calculated activationenergies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 404.2 Numbers of hydrogen bonds per water for the (6,6) nanotube . . . . . . . . . . 424.3 Radial pressures exerted on the (6,6) nanotube . . . . . . . . . . . . . . . . . . 434.4 Water self-diffusion coefficients and calculated activation energies . . . . . . . . 514.5 List of force field parameters of modified water models . . . . . . . . . . . . . . 525.1 Water/ion flow rates through a (9,9) nanotube and calculated activation energies 595.2 Ion transport efficiencies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 625.3 Water/ion entry rate into a (9,9) nanotube and calculated activation energies . 665.4 Water/ion diffusion coefficients and calculated activation energies . . . . . . . . 665.5 Percentage of ring-free water/numbers of hydrogen bonds per water for the(9,9) nanotube . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 686.1 Experimental data of real water and simulation results for different water models 76B.1 Effect of the cell dimension on water flow rates . . . . . . . . . . . . . . . . . . 99viiiList of TablesB.2 Effect of the thermostat employed on water flow rates . . . . . . . . . . . . . . 99B.3 Effect of the thickness of the force-exerted region on water flow rates . . . . . . 99B.4 Effect of the update frequency of the force exerted region on water flow rates . 99B.5 Effect of the external force on ion flow rates . . . . . . . . . . . . . . . . . . . . 100ixList of Figures2.1 Illustrations of water model geometries . . . . . . . . . . . . . . . . . . . . . . . 112.2 Illustration of chiral indices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122.3 Illustrations of a (6,6) carbon nanotube . . . . . . . . . . . . . . . . . . . . . . 122.4 Illustrations of the simulation cell . . . . . . . . . . . . . . . . . . . . . . . . . . 163.1 Water cumulative counts curves as functions of time through (8,8) and (9,9)nanotubes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 243.2 Snapshots of water configurations . . . . . . . . . . . . . . . . . . . . . . . . . . 253.3 Snapshot of the water structure in the (9,9) nanotube . . . . . . . . . . . . . . 263.4 Radial density profiles of water in (8,8) and (9,9) nanotubes . . . . . . . . . . . 293.5 Probability distributions for the number of water in (8,8) and (9,9) nanotubes . 303.6 Illustration of the cluster-by-cluster conduction mode . . . . . . . . . . . . . . . 333.7 Illustration of the diffusive conduction mode . . . . . . . . . . . . . . . . . . . . 344.1 Water cumulative counts curves as functions of time through the (6,6) nanotube 384.2 Snapshot of the water structure in the (6,6) nanotube . . . . . . . . . . . . . . 414.3 The dependence of ln(Rflow) on T−1 . . . . . . . . . . . . . . . . . . . . . . . . 454.4 The dependence of ln(Rentry) on T−1 . . . . . . . . . . . . . . . . . . . . . . . . 464.5 Rflow/Rentry dependence on T . . . . . . . . . . . . . . . . . . . . . . . . . . . . 474.6 Water potential energy profiles along the nanotube axis . . . . . . . . . . . . . 494.7 Water potential of mean force profiles along the nanotube axis . . . . . . . . . 494.8 The dependence of ln(Dwater) on T−1 . . . . . . . . . . . . . . . . . . . . . . . 504.9 RflowT/Dwater and RentryT/Dwater dependence on T . . . . . . . . . . . . . . . 53xList of Figures5.1 Water/ion cumulative counts curves as functions of time through the (9,9)nanotube . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 585.2 Ion flow rates as functions of the solution concentration . . . . . . . . . . . . . 605.3 The dependence of ln(Rflow,water) on T−1 . . . . . . . . . . . . . . . . . . . . . 645.4 The dependence of ln(Rflow,ion) on T−1 . . . . . . . . . . . . . . . . . . . . . . . 645.5 Rflow/Rentry dependence on T . . . . . . . . . . . . . . . . . . . . . . . . . . . . 675.6 Snapshots of ion environment in the (9,9) nanotube . . . . . . . . . . . . . . . 705.7 Sodium/Chloride ion potential of mean force profiles along the nanotube axis . 705.8 RflowT/Dwater and RentryT/Dwater dependence on T . . . . . . . . . . . . . . . 71xiList of Symbols and AbbreviationsA Area of MembraneD Diffusion CoefficientEa Activation EnergyNnt Number of Confined Water MoleculesP PressureR Gas ConstantRentry Entry RateRflow Flow RateRStokes Stokes RadiusT TemperatureU Potential EnergyVnt Volume of Carbon Nanotube Cavityd Bond Lengthdeff Effective Carbon Nanotube Diameterdnt Carbon Nanotube Diameterf Forcefex External ForcekB Boltzmann Constantkr Harmonic Force Constantn Number of Water Molecules being Applied External ForcexiiList of Symbols and Abbreviations(n,m) Chiral Indicesq Point Charger Coordinate Vectorr Radial Position inside Carbon Nanotubereff Effective Carbon Nanotube Radiust TimeuC Coulombic PotentialuLJ Lennard-Jones Potentialv Velocity Vectorw Potential of Mean Forceǫ0 Vacuum Permittivityε Lennard-Jones Energy Parameterθ Bond Angleη Fluid Viscosityξ Ion Transport Efficiency Parameterρ Particle Number Densityσ Lennard-Jones Length ParameterCNT(s) Carbon Nanotube(s)GROMACS Groningen Machine for Chemical SimulationsLINCS Linear Constraint SolverLJ Lennard-JonesMD Molecular DynamicsPME Particle Mesh EwaldPMF Potential of Mean ForceRC Reduced Charge ModelRP Reduced Lennard-Jones Parameter ModelxiiiList of Symbols and AbbreviationsSPC/E Extended Single Point Charge ModelTIP3P Three-site Transferable Intermolecular Potential ModelTIP4P/2005 Four-site Transferable Intermolecular Potential Model of 2005WHAM Weighted Histogram Analysis MethodxivAcknowledgementsFirst of all, I would like to express my sincere gratitude to my research supervisor, ProfessorGren Patey. I am deeply impressed by his knowledge as a scholar, his diligence as a learner,and his kindness as an instructor. He is like a blazing torch in the murky maze of science,illuminating me the path to the truth. I also appreciate his patient assistance over years inhelping me improving writing and presentation skills. I would not have accomplished thisdegree without his support.Besides my advisor, I am grateful to my fellow group members for their cooperation,encouragement, and friendship. My thanks go to Dr. Erin Lindenberg, Dr. Jingyi Yan, andDr. Sarah Overduin who gave me hands at the very beginning of my graduate study. I am alsograteful to my committee members and chemistry department staff for all kinds of assistance.Last but not the least, I would like to express profound gratitude to my parents for theiraltruistic support through my entire life. A special appreciation goes to Dr. Zhiwen Chenwho did me great favors during my study in Tianjin and Vancouver. I also would like tothank my friends all around the world for their warm companionship.xvChapter 1IntroductionThe transport of water and ions through nanoscopic channels occurs commonly in a vari-ety of physical systems. For example, water transport proteins, or aquaporins, are found inbiological membranes which can efficiently and exclusively conduct water molecules.1–4 Themolecular structures5–7 and water conduction mechanisms4, 8, 9 of aquaporins have been exten-sively illuminated by many modern techniques such as molecular dynamics (MD) simulationsand X-ray crystallography. Likewise, ion channels are those proteins which selectively con-vey ions across cell membranes.10–13 One major advantage of transport proteins is that theyconduct one specific particle or certain types of particles,2, 3, 10, 11, 14 avoiding the passage ofundesired molecules or ions. Also, growing research efforts demonstrate connections betweentransport proteins and various diseases.15–20 Therefore the study of transport proteins andtheir analogs will remain a topic of interest for decades.1.1 Molecule Transport through CNTsCarbon nanotubes (CNTs) are allotropes of carbon. Although CNTs are of much interest fortheir electrical,21 mechanical22 and thermal23 properties, in this dissertation we concentrate ontheir ability as conduction pores for water and ions. Conceptually, an open-ended, single-wallCNT can be constructed by “rolling up” a graphene sheet, which is simply one atomic layer ofgraphite. The synthesis methods of CNTs include arc discharge, laser ablation, chemical vapordeposition (CVD), et cetera.24 The lengths of synthesized nanotubes vary from the order ofmicrometers to millimeters, and the diameters of cylindrical single-wall nanotubes can beless than 1 nm, while the diameters of multiple-wall nanotubes are commonly several tens11.1. Molecule Transport through CNTsof nanometers.24 Resembling cavities inside transport proteins, the space surrounded by thecarbon scaffold allows the accommodation of molecules with diameters of several A˚ngstro¨mssuch as water,25–31 ethanol,32, 33 and benzene.34 Nanoscopic compounds, for example, single-stranded deoxyribonucleic acids (DNAs)35 and Iron(III) oxide nanoparticles,36 can be situatedwithin even larger CNTs.Although it is feasible to experimentally confirm the presence of water inside CNTs employ-ing various analysis instruments,25–31 it is still challenging to decipher complicated molecularstructures from spectra. Based on X-ray scattering data, radial density profiles of confinedwater have been calculated by Paineau et al.,30 indicating the presence of multiple water lay-ers. More recently Bernardina et al.,31 with the help of infrared spectroscopy, have suggestedthat confined water molecules provide dangling O−H bonds directed towards the carbon walls.Experimental scientists33, 37, 38 have also demonstrated that molecule transport through CNTsis enhanced compared to what one would expect based on classical fluid dynamics. Holt etal.37 have reported fast water and gas conduction rates through membranes which were madeof CNTs with diameters of less than 2 nm, and Qin et al.38 have discussed the dependence ofrate enhancement factors on the nanotube diameter. Moreover, excellent transport efficiencyacross nanotube membranes is not exclusive to water but common for various liquids, includ-ing ethanol and hexane.33 Efforts have been made on the application of the special dynamicalproperties of CNTs, for example in the realm of water desalination. A recent study39 hasshown that CNT-embedded polyamide membranes, which are widely used in reverse osmosis,can conduct water more rapidly than traditional, unmodified polymeric membranes.Furthermore, proof of ionic conductance through uncapped CNTs has been given in re-cent articles.35, 40–44 Ion transport can be detected by monitoring electric current tracechanges.40–43 Interestingly, the mobilities of alkali cations through isolated single-walledCNTs41, 42 and CNT-doped epoxy membranes44 have been proved very high, even exceed-ing their mobilities in bulk liquid. Unlike water, however, ion transport can be easily stoppedif nanotube terminals are chemically modified.4521.2. Computational Studies of Molecule Transport through CNTs1.2 Computational Studies of Molecule Transport throughCNTsUnlike experimental studies, which mostly investigate a large number of widely size-distributedand poorly characterized CNTs, computational studies focus more on a small number of chem-ically simple nanotubes. The structure of a simple CNT can be described by a set of integers(n,m) which are called chiral indices. The details of chiral indices are presented in Chapter2.2.Because of the nature of graphene,46 it is logical to expect that CNTs can not be wettedby water easily; yet simulation studies47–53 refute this assumption. An early study carriedout by Hummer et al.47 observed that water can automatically enter a short, narrow CNT(with an effective diameter of ∼0.5 nm). There is a debate about whether energy or entropydrives the entry of water. Kumar et al.48, 49 have claimed that water gains rotational entropyinside CNTs compensating for the increase in energy. Later Pascal et al.50 have concludedthat entropy is the primary driving force unless the diameters of CNTs are beyond a criticaldiameter of ∼1.1 nm. Studies51, 52 also have shown that entropy is the driving factor if CNTsare partially occupied, but energy is more favorable when CNTs are fully occupied. However,Waghe et al.53 have concluded that the filling is mainly induced by a decrease in energy,along with a minor increase in entropy.Unlike in bulk liquid, confined water can form special structures inside CNTs. Due tospatial constraints, in narrow CNTs, for example, a (6,6) nanotube, only a single molecularchain configuration is observed.47, 54–57 Water molecules are aligned one next to the other,connected by hydrogen bonds. Complex configurations develop as the CNT diameter in-creases. Using MD methods, Koga et al.58, 59 applied axial pressure on confined water at lowtemperature, and proposed a family of polygonal structures, where water molecules cluster aspolygonal structures with long-range order in the axial direction of intermediate-size CNTs(with effective diameters of ∼1.0 nm). Hydrogen bonds play a fundamental role in stabilizingring structures. Such quasi-one-dimensional structures were confirmed experimentally.29, 6031.2. Computational Studies of Molecule Transport through CNTsFor confined water in the liquid state, molecules can form similar polygonal (single layer)configurations but the axial periodicity vanishes.29, 55, 61–64 Other patterns have also beenreported,65, 66 where water molecules form helical chains along the CNT axis. Within largeCNTs (with effective diameters of ∼1.5 to 2.0 nm) confined water creates multilayer struc-tures,61, 67–69 but water inside even larger CNTs does not feel the nanotube wall, and thestructure resembles that of bulk liquid.61, 64, 68, 69On the macroscopic scale, pressure-driven liquid flow of an incompressible, Newtonianfluid inside a cylindrical pipe is often described as the Poiseuille flowu(r) =14η(dPdz)(r2eff − r2) , (1.1)where u is the flow velocity profile, η is the fluid viscosity, (dP/dz) is the pressure gradient, reffis the effective radius of pipe, and r is the radial position. The flow meets the no-slip conditionwhere the fluid has zero velocity at the fluid-solid boundary. Note that Equation (1.1) assumesthat the fluid is uniform and homogeneous. However, this assumption breaks down at thenanoscopic scale. Water molecules move cooperatively within narrow CNTs.47, 56 Hanasakiand Nakatani68 reported that water flow in intermediate-size CNTs has a flat velocity profilerather than a convex velocity profile, which is not compatible with Poiseuille flow. Theyfound that the longer lifetimes of hydrogen bonds between confined water leads to the fastconcerted motion of molecular clusters. Other papers55, 70 support this observation. The flowenhancement hints that the no-slip condition likely does not hold. MD studies71, 72 predict thatthe slip lengths of water decrease as CNTs become larger. The relative fast water conductionthrough CNTs is found to be closely related to the confined structure. Some articles63, 72 haveshown that increasing the curvature of CNTs can significantly reduce the friction at the fluid-solid interface. Joseph and Aluru73 concluded that the hydrophobic surfaces could influencewater molecule orientation in the proximity of CNTs and, therefore, enhance conduction rates.Several studies have discovered that water flow rates are sensitive to particular features ofCNT conduits, for example, the electrical charge.74–80 Most of the MD simulations mentioned41.2. Computational Studies of Molecule Transport through CNTsabove adopt the approximation that no carbon atoms have partial charges, yet this is not thecase.81, 82 Using ab initio MD simulations Won et al.74 and Sahu et al.75 investigated theeffect of partial charges on the dynamics of water. They concluded that partially chargedCNTs can induce a dipole moment, reorient water molecules, and make them more readilyaccessible to CNT orifices. Lu76 proposed that the dipole moment also prevents single-filemolecules from flipping, and, therefore, enhances water transport rates. Li et al.77 devised acontrollable nanoscopic fluid switch by adjusting the position of external charges. Other fac-tors that influence flow include conduit defects,83, 84 channel hydrophobicity,85–87 and chemicalmodification.75, 88, 89It is intuitively evident that the diameters of CNTs must exceed the diameters of ionsin order to allow ion entry and conduction. Using equilibrium MD simulations, Peter andHummer90 discovered that narrow CNTs block ions, but intermediate-size CNTs are Na+permeable. However, computer simulations91 suggested that CNTs can reject hydrated F−whose radius is greater than the ionic radius of F−. Theories about ion hydration/dehydrationare often mentioned to resolve such problems. In bulk solution, water molecules arrangeadjacent to ions, composing layers of hydration shells. The effective radius of a hydrationshell is determined by the ionic radius and the number of charges an ion bears. Studies92–94have demonstrated that the structure of hydration shells is strongly influenced by confinementin CNTs, therefore, to enter narrow CNTs some water molecules have to be removed from thehydration shell, giving rise to free energy penalties.90, 91, 95, 96 In contrast, within wide CNTs,ions are better stabilized by interacting with more neighboring water molecules.92 Generallyspeaking, it is less energetically expensive for ions to move through wide CNTs due to theminor alteration of hydration shells.Computer simulations have revealed that the CNT diameters strongly influence the dy-namics of ion transport through nanotubes. The scaffold of narrow CNTs can hinder thethree-dimensional ion diffusion.97–99 While enhanced ion mobilities in wide CNTs (with ef-fective diameters of ∼2.0 nm), compared to the bulk solution, have been reported,100 becauseof the recovery of bulk-like water structure but with less ion pairing. Moreover, other fac-51.3. Water Modelstors can also influence the ion dynamics, including charge distribution on CNTs,101–103 ionconcentration,104 and external electric fields.105CNTs have been suggested as a promising material for nano-filtration. As indicated above,narrow nanotubes can naturally exclude ions without obstructing water transport.90, 95, 106Hence nanotubes are promising candidates as nanoscopic filters for desalination of sea wa-ter. Although intermediate-size nanotubes permit ions to pass through, computational re-search88, 89 has demonstrated that chemical modification of such CNTs can optimize the effi-ciency of nanotubes as filtration systems. Also, it is possible to design nanoscopic devices toselect or detect specific ions. For example, the selection between Na+ and K+ using CNTsis a topic of interest. Investigations107–110 have revealed that nanotube preferences for differ-ent ions are primarily regulated by the particular hydration shell structures of Na+ and K+(including ion coordination number, the effective radius of hydration shell, and the intensityof which an ion holds surrounding waters) under specific conditions. Gong et al.111 pointedout that charge attraction/repulsion between modified CNTs and ions can improve selectionability as well.1.3 Water ModelsIn molecular simulations, a water model is a set of interaction parameters that describe thephysical properties of explicit water molecules in the aqueous phase. The first water modelwas proposed by Bernal and Fowler112 in 1933. The optimization of water models has beena work in progress for almost a century; several tens of water models have been proposed— the three-site transferable intermolecular potential model (TIP3P),113 the extended singlepoint charge model (SPC/E),114 the atomic multipole optimized energetics for biomolecularapplications model (AMOEBA),115 the four-site transferable intermolecular potential modelof 2005 (TIP4P/2005),116 and the six-site simple water model (SWM6),117 just to name afew.Unfortunately, there is not a single water model today which can perfectly describe all61.3. Water Modelsphysical properties of water in nature. In parameterizing water molecules, many factorsmust be considered all together, including quantum mechanics, molecular mechanics, andexperimental results.118 A water model is often detailed in three aspects. First, the numberof interaction sites. Considering Coulomb interactions, for the most widely used water modelssuch as TIP3P and TIP4P/2005, a point charge approximation is commonly applied, wheremolecules are described by several sites bearing either positive or negative partial charges.However, a limited number of sites poorly reproduce the actual charge distribution. An earlierstudy119 has demonstrated that a water molecule is better approximated by a combinationof point charges on nuclei and diffuse spherical charge densities around them. Second, themodel can be flexible or not. The use of rigid bonds where distances between interaction sitesare fixed to simplify computation complexity is often assumed, but in reality, chemical bondsare constantly vibrating. Models taking molecular flexibility into consideration can give moreaccurate results for some physical properties of water obtained in simulations, for example, theheat capacity at constant volume.120 To introduce flexibility in water models, the oscillationof bonds and angle can be described as simple harmonic potentials or more complicatedforms.121–123 Third, the model can be polarizable or not. In empirical models, it is oftenassumed that molecular polarizability is unaffected by the environment. For nonpolarizablewater models, partial charges are specifically parameterized to generate an appropriate dipolemoment114 similar to that of real water in liquid state. However, these models may not beaccurate in anisotropic environments, such as the proximity of a water-membrane interfaceor in the presence of an electric field. When inventing polarizable water models, the induced-dipole technique, the Drude-oscillator technique, and the fluctuating-charge technique arethree popular methods to explicitly take molecular polarizability into account.124An all inclusive, complicated water model can certainly enhance the accuracy of simulationresults, but can also considerably increase the computational cost. It is important to keepa balance between the efficiency of calculation and the reliability of conclusions. In thiscase, it is necessary to examine and compare different simplified water models and identifythose which will compromise least on simulation accuracy. Additionally, note that almost71.3. Water Modelsall water models are designed to reproduce the physical characteristics of bulk water, andit is likely that such properties can get skewed within confined regions, for instance withinCNT cavities. Therefore, deciding the best water model candidates for special anisotropicconditions is important.Investigations of the influence of different water models were pioneered by Alexiadis andKassinos.125–127 They systematically examined some water properties in single-wall CNTswith diameters ranging from 0.7 to 5.5 nm using MD simulation. Two families of three-site water models, TIP3P and SPC,128 were investigated, including both rigid and flexiblemodels. Although the confined water densities were similar within a (13,5) nanotube, they126suggested that the rigid SPC/E model produces an obvious pentagonal configuration, whilethe structure of the rigid TIP3P model is less ordered. They127 also discussed the effect offlexibility of water versions. By squeezing H−O−H bond angles, flexible water molecules arepacked more densely inside CNTs, and, therefore, greater average numbers of hydrogen bondswere observed in simulations, yet the variances of the numbers from rigid models were lessthan 10%. The flexible water models showed high similarities on selected physical properties,even though they are from two different families.Inspired by the early work, other studies have been conducted during the past five years.Nakamura and Ohno129 investigated the influence with more families of water models: three-site (TIP3P and SPC/E), four-site (the four-site transferable intermolecular potential model,TIP4P113), and five-site (the five-site transferable intermolecular potential model for usewith Ewald sum, TIP5P-E130). They investigated water configurations within (8,8) and (9,9)CNTs. It was shown that the water model employed can strongly influence the water struc-ture. TIP3P water is disordered, as noted in Reference 126; SPC/E water has structuresof stacked water clusters; TIP5P-E has complicated single-helix structures running throughCNTs. Moreover, the confined configurations of TIP4P water resemble those of SPC/E in the(8,8) CNT but are similar to TIP3P in the (9,9) CNT. Kumar et al.131 evaluated the effectsof model flexibility and polarizability. They discovered that the inclusion of polarizabilityimproved the accuracy of results for some dynamical properties (as compared to ab initio81.4. Outline of Dissertationsimulations), but that the inclusion of flexibility did not. They concluded that the classicSPC/E model appears to be the most excellent candidate for the study of confined watermolecules because it leads to sufficiently accurate results with lower computational cost.1.4 Outline of DissertationThe overall purpose of this dissertation is to study how and why different water modelsinfluence the simulated particle conduction through CNTs. It is interesting to see that bothwater and ion conduction rates are very sensitive to the water model employed in simulations.For water, we attribute such difference to two factors, bulk fluid diffusivity and molecularstructure within nanotubes. For ion transport, we provide evidence that the water structureinside an intermediate-size nanotube also has an impact on ion transport through CNTs.This dissertation is presented in six chapters. Chapter 1 overviews the background of wa-ter and ion conduction through CNTs and points out the necessity of studying water modelinfluences. Chapter 2 describes the molecular models employed, the simulation setup, andvarious analysis techniques. Chapters 3 and 4 thoroughly investigate, for both narrow andintermediate-size nanotubes, water conduction, confined water structure, and the related ther-modynamic and continuum hydrodynamic features that arise from different models. Chapter5 discusses the influence of water models on water and ion transport through an intermediate-size nanotube. Finally, Chapter 6 summarizes the conclusions and comments on the futureoutlook in this field.9Chapter 2Models and Methods2.1 Modelling of WaterIn this dissertation, we examine three widely used water models: the three-site transferable in-termolecular potential model (TIP3P),113 the extended single point charge model (SPC/E),114and the four-site transferable intermolecular potential model of 2005 (TIP4P/2005).116 Weselect the TIP3P and SPC/E models because they are among the most popular water modelsin use, and include the TIP4P/2005 model because it best reproduces the properties of ambi-ent bulk water.132 The first two models are three-site models, including one oxygen atom andtwo hydrogen atoms, while the third is a four-site model, with an additional virtual site atthe bisector of the H−O−H bond angle, representing electron distribution around the oxygenatom. The molecular geometries of the water models are displayed in Figure 2.1, and thegeometric parameters are summarized in Table 2.1. Note that all water models investigatedare rigid and nonpolarizable.Water model dO−H (nm) dO−M (nm) θH−O−H (degree)TIP4P/2005 0.0957 0.0155 104.5SPC/E 0.1000 n/a 109.5TIP3P 0.0957 n/a 104.5Table 2.1: List of geometric parameters of water models. d represents the distance betweensites, and θ represents the bond angle. The oxygen atom, hydrogen atoms, and the virtualsite are denoted as O, H, and M, respectively.102.2. Modelling of CNTsFigure 2.1: Illustrations of three-site (left) and four-site (right) water model geometries.Oxygen atoms are red, hydrogen atoms are gray, and the virtual midpoint site is pink.2.2 Modelling of CNTsThe structure of a CNT can be described by a set of integers (n,m) which are called chiralindices. The set defines how the nanotube is “rolled up” from a graphene sheet (see Figure2.2). If n = m the CNTs are called armchair nanotubes, and if m = 0 they are called zigzagnanotubes. Otherwise, CNTs are termed chiral. Figure 2.3 illustrates a pristine (6,6) CNT.It is possible to calculate the diameter of a nanotube, dnt, from its chiral indices (n,m)dnt =aπ√3(n2 + nm+m2) , (2.1)where a is 0.142 nm, the length of a carbon-carbon bond. An effective diameter, deff , can alsobe calculated if we take the van der Waals radius of a carbon atom, rcarbon, into considerationdeff = dnt − 2rcarbon , (2.2)where dnt is from Equation (2.1) and rcarbon is 0.17 nm.95 Note that in the current investigationwe assume that CNTs are rigid with no termination, which is an approximation regularly usedin CNT studies.47, 77, 88, 95, 97, 107, 129 MD research articles investigating the effects of nanotube112.2. Modelling of CNTsFigure 2.2: Illustration of chiral indices. a1 and a2 are unit vectors of the graphene sheet. An(n,m) CNT is conceptually constructed by “rolling up” the graphene plane along the vectorCh= na1+ma2; T denotes the CNT axis. This figure is from Wikipedia with permission.Figure 2.3: Illustrations of a side view (left) and a top view (right) of a pristine (6,6) CNT.Carbon-carbon bonds are depicted as cyan cylindrical segments.122.3. Interaction Potentialsflexibility92, 126, 127, 133 and termination75, 89 are available.In this dissertation, the terms CNT(s) and nanotube(s) are used interchangeably.2.3 Interaction PotentialsIn our models, all site-site interactions involve Lennard-Jones (LJ) and Coulombic terms. TheLJ potential is a model to approximate the van der Waals interaction between a pair of atoms.uLJ,ij, the LJ potential between atom i and j, has a simple expressionuLJ,ij = 4εij[(σij|rij |)12−(σij|rij |)6], (2.3)where εij is the depth of the potential well, σij is the distance at which uLJ,ij is exactly zero,and |rij | is the distance between atoms i and j. The pair LJ parameters (εij and σij) arecomputed by the Lorentz-Berthelot combination rule from atomic force field parametersεij =√εiεj , (2.4)σij =σi + σj2. (2.5)In practice, given that uLJ,ij converges to zero quickly, it is reasonable to presume that thelong-range uLJ,ij is negligible. The LJ cut-off radius in this dissertation is 1.2 nm, and thetotal LJ contribution to the potential energy will be simply the sum of all pairwise interactionswhose |rij | is within this arbitrary cut-off radius. The pair electrostatic potential uC,ij of sitesi and j who both bear charges is given by Coulomb’s LawuC,ij =qiqj4πǫ0|rij | , (2.6)where qi and qj are charges on interaction sites i and j, and ǫ0 is the vacuum permittivity.However, unlike uLJ,ij, uC,ij decays to zero very slowly, which means the long-range uC,ij still132.4. Simulation Detailssignificantly contributes to the total potential. Given that the periodic boundary conditionsare used in MD simulations, as a result, we have to consider interactions of a given chargewith not only other charges in the same simulation cell but also with all images (includingitself) in the periodic cells. In practice, the particle mesh Ewald (PME) method134, 135 is usedto calculate the total electrostatic interactions. In this method, the slowly converging sum canbe obtained from two quickly converging terms, one in real space and the other in reciprocalspace. The real space PME electrostatics is truncated at 1.2 nm. The reciprocal space PMEsummation uses more than 85 wave vectors in each direction.All atomic force field parameters are summarized in Table 2.2. Parameters for a car-bon atom are from the AMBER03 force field.136 Parameters for different water models arefrom References 113, 114, and 116. Parameters for ions are those proposed by Joung andCheatham,137 which are specifically optimized for different water models. Because optimizedion parameters for the TIP4P/2005 model are not available, we use those for the TIP4P-Ewmodel instead. Furthermore, two arbitrarily modified TIP4P/2005 water models are designedin this dissertation so as to reveal the mechanism by which water transport through nanotubesdepends on the water model employed. The revised force field parameters are described inChapter 4.4.2.4 Simulation DetailsThe simulation cell is illustrated in Figure 2.4. The cell is a rectangular hexahedron withCartesian dimensions (x, y, z) of (5.116 nm, 5.168 nm, 7.569 nm) unless otherwise stated.Two rigid graphene sheets with appropriate openings are located in xy planes, and an open-ended CNT, whose symmetry axis is along the z direction, is embedded within the graphenesheets. In this dissertation, we examine three sizes of armchair CNT: (6,6), (8,8), and (9,9).The length of the nanotube is fixed at 3.561 nm. From Equation (2.2) the effective diametersof (6,6), (8,8) and (9,9) CNTs are 0.474, 0.745, and 0.881 nm, respectively. The skeletonof carbon atoms can be regarded as a permeable membrane, separating the cell into two142.4. Simulation DetailsAtom σ (nm) ε (kJ mol−1) q (e)CNTC 0.3400 0.3598 0.0000Water — TIP4P/2005O 0.3159 0.7749 0.0000H 0.0000 0.0000 +0.5564M 0.0000 0.0000 −1.1128Water — SPC/EO 0.3166 0.6502 −0.8476H 0.0000 0.0000 +0.4238Water — TIP3PO 0.3151 0.6364 −0.8340H 0.0000 0.0000 +0.4170Ion — TIP4P-EwNa+ 0.2184 0.7050 +1.0000Cl− 0.4918 0.0488 −1.0000Ion — SPC/ENa+ 0.2160 1.4761 +1.0000Cl− 0.4830 0.0535 −1.0000Ion — TIP3PNa+ 0.2439 0.3660 +1.0000Cl− 0.4478 0.1490 −1.0000Table 2.2: A list of the force field parameters used in the present dissertation.152.4. Simulation Detailszx / yyxFigure 2.4: Illustrations of a side view (top) and a top view (bottom) of the simulationcell, representing a permeable membrane which contains a (6,6) CNT connecting two waterreservoirs. Carbon atoms are cyan, oxygen atoms are red, and hydrogen atoms are gray. The(x, y, z) dimensions of the rectangular hexahedral cell are (5.116 nm, 5.168 nm, 7.569 nm).The region highlighted by blue shadow denotes the region where the external force is appliedto water oxygen atoms.162.5. Application of Pressure Differencereservoirs. During simulations, all carbon atoms are frozen at the initial positions.The reservoirs are filled with either water or NaCl solutions of varying concentration.The density of water was set at approximately 1 g cm−3, taking the van der Waals radiusof carbon atom into account when calculating the fluid volume. The systems included 3266,3299, and 3322 water molecules for (6,6), (8,8), and (9,9) CNTs, respectively. NaCl solutionswere considered only for the (9,9) nanotube. Initially, identical systems were replicated aspreviously described for the (9,9) nanotube; then water molecules were randomly substitutedwith either Na+ or Cl− ions. The simulations considered contained 15, 60, and 166 ion pairs,giving NaCl solutions of approximately 0.25, 1, and 2.8 mol L−1, respectively.We employed the GROMACS package138 (version 4.5.5, double precision) for MD sim-ulations. Periodic boundary conditions were implemented in all three dimension. All MDsimulations were carried out in the canonical (NV T ) ensemble using a time step of 2 fs. Mostsimulations were 10 ns in duration, with first 1 ns being sufficient to establish steady stateflow, and data were collected from the last 9 ns. For each system, at least five individualsimulations with different starting configurations were conducted, and these were used to cal-culate averages and estimated standard deviations. The system temperature was maintainedusing the velocity rescaling algorithm of Bussi et al.139 A brief description of the thermostatis given in Appendix A.3. We investigated fluid structural and dynamical properties at fourdifferent temperatures: 260, 280, 300, and 320 K.2.5 Application of Pressure DifferenceTo maintain a stable fluid flow through a CNT channel, a hydrostatic pressure difference,∆P , was applied. In this dissertation we used the algorithm proposed by Zhu et al.140 In thismethod, a constant force along the z axis is applied to a subset of water molecules. Then thepressure difference in z direction is determined by∆P =nfexA, (2.7)172.6. Data Analysiswhere n is the total number of molecules in the subset, fex is the constant external force, andA is the area of the membrane.The original algorithm proposed by Zhu et al.141 uses the same expression as Equation(2.7) to calculate the pressure difference but applies the force to all water molecules in bothreservoirs, including those adjacent to the orifice of the CNT conduit. However, this setupmay artificially accelerate the rate of water permeating the channel. This drawback can beovercome by applying the external force only to water molecules located at the top and bottomof the periodic simulation cell (see Figure 2.4). Suk and Aluru142 demonstrated that therevised algorithm could generate the correct pressure difference across the membrane. Onlythe oxygen atoms of water were subjected to the external force, avoiding instigating additionalmolecular rotation. In investigating the NaCl solutions, no force was applied to the ions, whichprevents the artifact of unequal concentrations between reservoirs. The regions where the forcewas exerted were 0.2 nm thick, including about 360 water molecules in the subset. The area ofthe membrane was 26.44 nm2 (5.116 nm × 5.168 nm). In implementing this algorithm on theGROMACS platform efficiently, the subset list was updated every 10 ps, rather than at everytime step. We verified that within reason altering the thickness of the regions or the periodof subset update has no significant qualitative effect, and some supplemented information ispresented in Appendix B.3. In our simulations, the target pressure difference was set at 220MPa, except when otherwise specified. Note that this pressure difference is too high to bephysically realistic, but a large value is necessary to ensure efficient samples on simulationtime scales.Other algorithms for providing a pressure gradient include a reflecting particle mem-brane,63 a fluidized piston model,68 a movable piston model.1432.6 Data AnalysisDuring the production phase of simulations, trajectories and velocities of molecules wereperiodically recorded. By examining and interpreting these results using programs and scripts182.6. Data Analysiswe obtain useful structural and dynamical information.2.6.1 Flow Rates and Entry RatesThe flow rate, Rflow, is defined as the number of particles that pass through the nanotube(entering from one end and leaving the other) per unit time. This definition is widely used insimilar investigations.73, 74, 85, 88, 95 Under a pressure difference of hundreds of MPa, transportevents are only observed in the direction of the external force.The entry rate, Rentry, is defined as the number of particles that enter from the feedreservoir into the nanotube interior per unit time. By examining the entry rate, we learn howfrequently a particle enters the nanotube orifice from the bulk liquid. For calculation purposes,a particle is considered to have entered the nanotube if its center passes the nanotube orificeshown in Figure 2.4. Note that for water the molecular center is taken to be the oxygen atom.Once entered, a particle can continue to traverse the nanotube or be rejected back into thebulk. To avoid overcounting entries due to the same particle “vibrating” back and forth acrossthe plane defining entry, we count only the first entry of any given molecule in the calculationof the entry rate.2.6.2 Characterization of Confined ConfigurationsThe formation of hydrogen bonds are ubiquitous between water molecules, and is a straight-forward property to investigate. We employ the hydrogen bond criteria proposed by Luzar andChandler.144 The criteria for a hydrogen bond are twofold: the angle of Oacceptor−Odoner−Hmust be less than 30 degrees and the distance between Oacceptor and Odoner less than 0.36 nm.Average numbers of hydrogen bonds for confined molecules were calculated during the dataanalysis.As noted in Chapter 1, different polygonal water clusters are commonly observed withinintermediate-size nanotubes. We identify just stacked ring structures in our simulations, butno long-range spiral configurations are found inside either (8,8) or (9,9) nanotubes. Here wedescribe the empirical method used to analyze ring clusters in detail. (1) Sort all oxygen192.6. Data Analysisatoms in ascending z coordinate in the interior of the nanotube. (2) Starting from the top ofthe list, define cluster boundaries between atoms whose z coordinate gap are greater than 0.1nm. (3) Calculate the z coordinate distance from the top to the bottom atom of each cluster.In our analysis, water molecules are considered to be part of the same ring configuration ifthe differences in the z coordinate of their oxygen atoms are all less than 0.15 nm, and thosemolecules are designated as ring-bound. We assign clusters to different categories, square,pentagonal, and hexagonal corresponding to those having 4, 5, and 6 molecules, respectively.All remaining water molecules are labeled as ring-free. Note that the distance criteria (0.1and 0.15 nm) employed in the procedure as mentioned earlier are arbitrary, but they aresufficient to ensure that molecules are arranged as rings rather than small coils or randomconfigurations.2.6.3 Diffusion CoefficientTo better understand the mechanism of model dependence, we investigate bulk self-diffusioncoefficients for water and ions (Dwater and Dion). The MD simulations used for diffusioncoefficient calculations were conducted under equilibrium conditions in a fixed cubic cell oflength 3.000 nm. For pure water, there were 903 water molecules within the cell, such thatthe density of bulk liquid was 1.00 g cm−3. For NaCl solutions, random water molecules werereplaced by 4, 16, and 45 pairs of Na+ and Cl− ions, giving 0.25, 1, and 2.8 mol L−1 solutions,respectively. Each simulation was run for 100 ns at a constant temperature, the first 5 ns forequilibrium, and the last 95 ns for data collection. The diffusion coefficients, D, are obtainedusing the Einstein relation via mean square displacements in the usual manner52, 55D = limt→∞〈|r(t)− r(0)|2〉6t, (2.8)where |r(t) − r(0)| is the three-dimensional particle displacement during a time period of t,and the angular brackets indicate an average over all particles.202.6. Data Analysis2.6.4 Potential of Mean ForceThe calculation of free energy differences is essential in computational science. The potentialof mean force (PMF) is helpful in calculating free energy changes along particular reactioncoordinates. The basic concepts of PMF are presented in Appendix A.4. In this dissertation,we investigate the free energy profiles of water/ions along the symmetry axis of a CNT.To obtain the PMF profiles, w(z), for water, we use the average densities from equilibriumsimulations95w(z) = −kBT ln[ρ(z)ρ0], (2.9)where ρ(z) is the number density profile as a function of z coordinate along the axis, and ρ0is the density of bulk water.Free energy profiles for ions are calculated using umbrella sampling.145, 146 A bias potentialfunction, Ubiasi (z), of the harmonic form, was applied to a selected ionUbiasi (z) =12kr(z − zi)2 , (2.10)where kr is the force constant taken to be 1000 kJ mol−1 nm−2, and zi is the target zcoordinates along the central axis. The target positions were moved from 0.5 nm to 3.5 nmwith an interval of 0.1 nm. Each window simulation had a duration of 10 ns and a time stepof 2 ps. The weighted histogram analysis method (WHAM) algorithm,147 which calculatesPMF profiles, was achieved using the GROMACS built in program g wham.21Chapter 3Simulated Water Conduction Ratesthrough Intermediate-size CNTs ∗3.1 OverviewA previous MD study by Nakamura and Ohno129 investigated water structure inside (8,8) and(9,9) CNTs and discovered convincing structural differences amongst water models. However,their simulations and model comparisons apply only to equilibrium conditions, and not tothe nonequilibrium cases. Motivated by this work, we investigate whether confined waterstill exhibit model dissimilarities when undergoing pressure-driven flow. We find that theflow rate is significantly influenced by the water model employed. We trace the dynamicaldiscrepancies to the structure taken on by the water molecules inside nanotubes. Two distinctconduction modes are proposed to explain why flow rates through intermediate-size nanotubesare remarkably different for various water models.3.2 Water Flow Rates through Intermediate-size CNTsIn Figure 3.1 we plot typical cumulative counts of water molecules passing through both (8,8)and (9,9) nanotubes at 300 K as functions of time. At the pressure difference used in oursimulations, water only traverses the nanotube in the same direction as the external force.It is interesting to observe that, starting at 1 ns, the counts are essentially linear in timeregardless of the size of nanotube or the water model employed in the simulations. Figure 3.1∗A version of this chapter has been published. L. Liu and G. N. Patey, J. Chem. Phys., 141, 18C518(2014).148223.3. Water Structure inside Intermediate-size CNTsshows that a 1 ns period is sufficiently long to establish steady state flows. The average flowrate Rflow is commonly used to describe the conduction ability of CNTs. The different slopesrepresent the flow rates and are obtained by linearly fitting the cumulative counts curves.From Figure 3.1, it is also obvious that the flow rates are not identical for different watermodels.Average flow rates, which means the average number of water molecules passing throughthe nanotube during a time period of 1 ns, with estimated standard deviations at 300 K, forall three models and different nanotubes are tabulated in Table 3.1. The results show that thewater conduction rate is strongly dependent on the type of model employed. The TIP4P/2005model has the slowest flow rates while the TIP3P model has the fastest; flow rates for SPC/Eare intermediate, but closer to those of TIP4P/2005. The rate ratio between the slowest andfastest models is approximately five.3.3 Water Structure inside Intermediate-size CNTsBefore discussing our results, we briefly review the work of Nakamura and Ohno.129 Underequilibrium conditions, these authors studied the effect of different water models on molecularstructures within armchair CNTs considering a total of four water models: TIP3P, SPC/E,TIP4P, and TIP5P-E. Their MD simulations were conducted at 280 and 300 K. From densityprofiles, they found that the difference between models is negligible inside narrow nanotubes,but becomes significant inside intermediate-size CNTs. In particular, for (8,8) and (9,9)nanotubes, confined TIP3P water has a lower density, and less ordered structure comparedWater model Rflow,(8,8) (ns−1) Rflow,(9,9) (ns−1)TIP4P/2005 51 (4) 107 (6)SPC/E 111 (4) 134 (14)TIP3P 287 (5) 492 (15)Table 3.1: Average flow rates for different water models in (8,8) and (9,9) nanotubes. Thesimulations were conducted at 300 K. The numbers in parenthesis are estimated standarddeviations.233.3. Water Structure inside Intermediate-size CNTs0 2 4 6 8 10Time (ps)050010001500200025003000Cumulative CountsTIP4P/2005SPC/ETIP3P(a)0 2 4 6 8 10Time (ps)010002000300040005000Cumulative CountsTIP4P/2005SPC/ETIP3P(b)Figure 3.1: Cumulative counts curves showing the number of water molecules that passthrough (8,8) [panel (a)] and (9,9) [panel (b)] nanotubes as functions of time. Note thatcounts during first 1 ns are not displayed because this is the equilibration period. The red,blue, and green curves are for the TIP4P/2005, SPC/E, and TIP3P models, respectively.243.3. Water Structure inside Intermediate-size CNTsFigure 3.2: Examples of local water configurations inside nanotubes: a square configuration ofTIP4P/2005 water in the (8,8) nanotube (top left), a pentagonal configuration of TIP4P/2005water in the (9,9) nanotube (top right), a hexagonal configuration of SPC/E water in the (9,9)nanotube (bottom left), and a less ordered TIP3P water cluster (bottom right). Carbon atomsare cyan, oxygen atoms are red, and hydrogen atoms are gray. The black dotted lines indicatehydrogen bonds.to the other water models considered. The SPC/E and TIP4P models result in stacked waterrings structures, while TIP5P-E water molecules form a spiral structure running throughthe nanotube. They investigated confined water orientation as well. For the SPC/E andTIP4P models, they observed antiferroelectric arrangements with the net dipole moment alongnanotube axis fluctuating about zero. For the TIP5P-E model, they observed a ferroelectricarrangement with a non-zero dipole moment in the axial direction. Finally, they also predictedthat for SPC/E and TIP4P water, ring structures would have to enter/exit nanotubes as awhole. Because of the structural diversity under equilibrium conditions, we were curious asto whether or not such diversity remains under nonequilibrium circumstances.In Figure 3.2 we illustrate examples of local water structures inside nanotubes, including253.3. Water Structure inside Intermediate-size CNTsFigure 3.3: A configurational snapshot of TIP4P/2005 water in the (9,9) nanotube. Carbonatoms are cyan, stacked pentagonal ring structures are repeating green, orange and yellow,ring-free water molecules are magenta. Note that part of the carbon nanotube is not displayedfor visual clarity.different ring configurations and less ordered clusters. Water molecules which form a ringconfiguration (“ring-bound” as defined in Chapter 2.6.2) associate with hydrogen bonds. Ow-ing to spatial constraints, there are very few pentagons, and no hexagons found in the (8,8)nanotube, whereas, all three ring structures do exist in the (9,9) nanotube. Water moleculesfrom less ordered clusters (“ring-free” as defined in Chapter 2.6.2) appear as more randomlypositioned and oriented. A configurational snapshot of confined water in a (9,9) nanotube isshown in Figure 3.3. Stacked ring structures are highlighted using repeating green, orange,and yellow color schemes. The disordered molecules are colored magenta. It is also interestingto note that the stacked ring clusters of liquid water, to some extent resemble the ice struc-tures reported in CNTs at lower temperatures.58, 59 Nevertheless, axial long-range order, asfound for ice, does not occur because of the presence of ring-free molecules together with therandom molecule orientation of individual clusters. Note that molecules can and do switchbetween ring-bound and ring-free states as the simulation evolves.To gain insight into how the water structure varies for different models in (8,8) and(9,9) CNTs, we quantitatively inspected all systems considered for ring-bound and ring-freemolecules; to be more specific, we divided the ring-bound molecules into those that are mem-bers of square, pentagonal, or hexagonal rings. The method employed is described in Chapter2.6.2. The results are summarized for both nanotubes in Table 3.2. It is apparent that263.3. Water Structure inside Intermediate-size CNTsWater model Square (%) Pentagon (%) Hexagon (%) Ring-free (%)(8,8) CNTTIP4P/2005 96.7 0.2 0 3.1SPC/E 74.2 0.6 0 25.2TIP3P 28.7 0.3 0 71.0(9,9) CNTTIP4P/2005 3.3 49.2 34.9 12.6SPC/E 2.8 12.3 72.3 12.6TIP3P 12.9 20.0 3.0 64.1Table 3.2: The average percentage of water molecules in different structural states inside (8,8)and (9,9) nanotubes.the fractions of ring-bound and ring-free molecules differ considerably from one model to an-other. Amongst the three models considered, TIP3P water has the largest number of ring-freemolecules, with 71.0% in the (8,8) and 64.1% in the (9,9) nanotube, respectively. Our resultsquantitatively demonstrate less ordered structures for the confined TIP3P model, consistentwith the conclusion of Nakamura and Ohno.129 For the TIP4P/2005 and SPC/E models, localring structures are confirmed by relatively large percentages of ring-bound molecules. Withinthe (8,8) nanotube the square configuration is predominant. Interestingly, within the (9,9)nanotube, TIP4P/2005 favors the pentagonal configuration over the hexagonal (49.2% versus34.9%), while SPC/E is opposite (12.3% versus 72.3%). The configurational preference of theSPC/E model we obtain again coincides with the observations of Nakamura and Ohno.129It is evident that organized stacked ring structures will yield a notably sharp water radialdensity distribution profile, ρ(r). Density profiles at 300 K are plotted in Figure 3.4. We seethat the TIP3P model has the broadest distributions, especially in the (9,9) nanotube, wherethe non-zero ρ(r) near the central axis clearly reflects the existence of ring-free molecules.For the (8,8) nanotube, TIP4P/2005 water has a taller and narrower distribution curve thanSPC/E, which can be explained by its smaller percentage of ring-free molecules (3.1% versus25.2%). For the (9,9) nanotube, the TIP4P/2005 model favors smaller pentagonal rings ratherthan hexagon rings in contrast with SPC/E. Therefore, it is reasonable that the distribution273.3. Water Structure inside Intermediate-size CNTscurve of TIP4P/2005 shifts inward, and has a lower peak than SPC/E.We also considered the number of water molecules, Nnt, confined inside CNTs. Modeldifferences are also evident in probability distribution functions, P (Nnt), as shown in Figure3.5. First, we considered the (8,8) nanotube. For the SPC/E and TIP3P models, we note thatthe P (Nnt) distributions are unimodal, whereas for the TIP4P/2005 model the distributionis bimodal, with quite sharp peaks at 48 and 52 molecules. As 96.7% of the TIP4P/2005water molecules within the nanotube are in square ring configurations, the two peaks likelycorrespond to 12 and 13 square rings, respectively, stacked along the nanotube axis. Becausethe SPC/E model has the largest fraction (72.3%) of hexagonal rings in the (9,9) nanotube,we believe that the peak in P (Nnt) at 72 molecules comes from the contribution of 12 hexago-nal rings. The TIP4P/2005 model has significant numbers of both pentagonal and hexagonalrings (49.2% and 34.9%), as well as ring-free molecules, hence the number distribution forTIP4P/2005 is broad with a major peak at 64 molecules, implying mixed structures, and aminor peak at 72 molecules, which again suggests 12 hexagons. Moreover, inside both nan-otubes the TIP3P model generally contains the fewest numbers of molecules, which suggeststhat TIP3P water is the most loosely packed within the nanotubes.On account of the distinctive structures observed for the different models, it is interestingto investigate the average number of hydrogen bonds per water molecule (hydrogen bondnumber) within the nanotubes. Analogous to ice structure inside nanotubes,58, 59 if stackedring configurations are prevailing, we would expect both intra-ring hydrogen bonds gathermolecules together, and inter-ring hydrogen bonds connect neighboring rings. In contrast, ifthe water molecules are more loosely packed, we would predict on average fewer hydrogenbonds per molecule. The method of determining hydrogen bonds is described in Chapter2.6.2. The results are listed in Table 3.3 and are in agreement with our expectations. In bothnanotubes, TIP3P water has the smallest hydrogen bond numbers, in accord with its havingthe most substantial fraction of ring-free molecules. TIP4P/2005 water has on average morehydrogen bonds within the (8,8) nanotube than SPC/E (3.28 versus 3.02 per molecule), whichcan be explained by its higher fraction of ring-bound molecules (96.9% versus 74.8%). In the283.3. Water Structure inside Intermediate-size CNTs0 0.1 0.2 0.3 0.4r (nm)00.040.080.12ρ(r)TIP4P/2005SPC/ETIP3P(a)0 0.1 0.2 0.3 0.4 0.5r (nm)00.040.080.12ρ(r)TIP4P/2005SPC/ETIP3P(b)Figure 3.4: Water radial density profiles measured from the nanotube symmetry axes (r = 0)at 300 K inside (8,8) [panel (a)] and (9,9) [panel (b)] nanotubes. The water models areindicated in the legends.293.3. Water Structure inside Intermediate-size CNTs40 50Nnt00.10.20.30.4P(N nt)TIP4P/2005SPC/ETIP3P(a)60 70 80Nnt00.050.10.150.20.25P(N nt)TIP4P/2005SPC/ETIP3P(b)Figure 3.5: Probability distributions for the number of confined water molecules within (8,8)[panel (a)] and (9,9) [panel (b)] nanotubes. The water models are indicated in the legends.303.4. Water Transport MechanismsWater model (8,8) CNT (9,9) CNTTIP4P/2005 3.28 3.22SPC/E 3.02 3.28TIP3P 2.66 2.75Table 3.3: Average numbers of hydrogen bonds per water molecule (hydrogen bond numbers)inside the CNTs.(9,9) nanotube both TIP4P/2005 and SPC/E have the same fraction of ring-free molecules(12.6%), and therefore similar hydrogen bond numbers.3.4 Water Transport MechanismsThe results shown in Tables 3.1 and 3.2 indicate a clear correlation between the flow rateand the percent of ring-free molecules. In both nanotubes, TIP3P has the largest faction ofring-free molecules (71.0% and 64.1%) as well as the highest flow rates (287 and 492 ns−1).TIP4P/2005 water has fewer ring-free molecules (3.1%) for the (8,8) nanotube than SPC/E(25.2%) and a slower flow rate (51 versus 111 ns−1). For the (9,9) nanotube TIP4P/2005 andSPC/E have similar flow rates (107 versus 134 ns−1) consistent with identical fractions of ring-free molecules (12.6%). We do note that although the fraction of ring-bound molecules is thesame for both these models, the ring configuration distributions differ, indicating differencesin the structure of the confined water. However, these structural differences appear to have,if any, a minor effect on the flow rates.By carefully investigating simulation trajectories of water flow through nanotubes, we canindentify two conduction modes. The first mode is a “cluster-by-cluster” mode, which is bestrepresented by the TIP4P/2005 model. In this mode, water molecules aggregate into clusters(square, pentagonal or hexagonal rings) which move together through the nanotube, onecluster following the other. The cluster-by-cluster mode is depicted by simulation snapshotsfor TIP4P/2005 water in the (8,8) nanotube in Figure 3.6. From Figure 3.6 we can clearlyidentify multiple square ring configurations as discussed in Chapter 3.3. We point out that, inthis particular example, the highlighted rings existing at t = 0 ps remain intact as they travel313.4. Water Transport Mechanismsthrough the nanotube. We observe that the ring-bound and ring-free molecules may switchfrom one type to the other, therefore, rings do not always remain intact, and sometimesclusters do break apart as they flow through the nanotube. Still, this conduction mode isvalidated by many simulation trajectories, which show a significant amount of cluster motion.The second mode is a “diffusive” mode, which is the primary mode for the TIP3P model.In this case, there are fewer ring structures (Table 3.2), and those that do occur are fragileand tend to break up during transport, partially because of fewer hydrogen bonds (Table3.3). Therefore, water molecules pass through the nanotube more independently and freely,rather than as part of a cluster or group. The diffusive mode is illustrated in Figures 3.7with snapshots from simulations of TIP3P water in the (8,8) nanotube. Here we note thatthe apparent three square rings present at t = 0 completely break apart as they traversethe nanotube. We do not observe TIP3P clusters pass through the entire nanotube intact,suggesting that distinct cluster motion for TIP3P is very rare.Based on this analysis, we conclude that the formation of relatively stable stacked ringstructures for the TIP4P/2005 and SPC/E models leads to flow rates that are much slowerthan those observed for TIP3P water. The cluster-by-cluster mode requires concerted motionof molecules from several clusters, while the diffusive mode allows individual molecules tomove forward without substantial restriction. The different time scales also illustrate thedifferent flow rates of the two conduction modes in Figures 3.6 and 3.7. Of course, bothconduction modes can and do happen to some extent for all three models, but TIP4P/2005and SPC/E favor the slower cluster-by-cluster mode, where molecules are likely to assembleinto ring configurations, and TIP3P favors the faster diffusive mode, where there is a higherpercentage of ring-free molecules.It is perhaps worth mentioning that density differences alone can not account for theobserved differences in flow rates. This factor is most evident in the (9,9) case, where themost probable nanotube densities are similar for TIP4P/2005 and TIP3P models, as shownin Figure 3.5, which have very different flow rates, and significantly different for TIP4P/2005and SPC/E models, which have similar flow rates.323.4. Water Transport Mechanisms(a) t = 0 ps(b) t = 50 ps(c) t = 100 psFigure 3.6: Illustrations of the cluster-by-cluster conduction mode for TIP4P/2005 water inthe (8,8) nanotube. Carbon atoms are cyan, oxygen atoms are red, and hydrogen atoms aregray, except for three square ring configurations highlighted in yellow, orange, and green. Weremark that the highlighted ring structures pass intact across the nanotube. Note that partof the carbon nanotube is not displayed for visual clarity.333.4. Water Transport Mechanisms(a) t = 0 ps(b) t = 5 ps(c) t = 10 psFigure 3.7: Illustrations of the diffusive conduction mode for TIP3P water in the (8,8) nan-otube. Carbon atoms are cyan, oxygen atoms are red, and hydrogen atoms are gray, exceptfor three square ring configurations highlighted in yellow, orange, and green. We remark thatthe highlighted ring structures existing at t = 0 ps break up completely as the water moleculespass through the nanotube. Note that part of the carbon nanotube is not displayed for visualclarity.343.5. Summary3.5 SummaryWe have investigated pressure-driven water transport through (8,8) and (9,9) CNTs using MDsimulations. Among the three water models examined, TIP3P water has the highest transportrates in both nanotubes, which are about five times faster than those of the slowest model,TIP4P/2005. Next, we trace the strikingly different flow rates to the varying amounts of wateroccurring as ring-bound and ring-free molecules. The TIP4P/2005 model tends to constructstacked ring configurations, giving rise to a cluster-by-cluster conduction mode with many ringclusters traveling as single units through the nanotube. This is also true for the SPC/E model,but to a lesser extent. In contrast, the TIP3P model favors a diffusive conduction mode, wherering structures occur less frequently and are likely to break apart as water passes throughthe nanotube. The diffusive mode is faster than the cluster-by-cluster mode because watermolecules move as individual particles in the diffusive mode rather than as parts of largermolecular groups. Our simulations also support the conjecture of Nakamura and Ohno129who suggested that water rings could be regarded as a whole in terms of hydrodynamics.35Chapter 4How Water Conduction Ratesthrough CNTs Is Related to BulkProperties for Different WaterModels∗4.1 OverviewIt is interesting to notice that water conduction rates through intermediate-size CNTs arestrongly influenced by the water model employed.148 In Chapter 3 we traced such differ-ences to confined water structures and proposed two distinct transport modes to explain ourobservations. It is important to emphasize that the results for intermediate-size nanotubesdiscussed above can not be presumptuously extended to smaller nanotubes. For example, ina (6,6) nanotube, only simple single-file conduction with no complex ring configurations isobserved.47, 70, 95 However, we show that the flow rate discrepancies for different water modelsremain in the (6,6) case. By carrying out simulations at various temperatures and carefullyanalyzing the water dynamics, we show that the discrepancies, in fact, reflect the differentmobilities of the bulk liquids. Also, our results are consistent with continuum hydrodynamicanalysis. We revisit water flow rates through a (9,9) nanotube. Our work illuminates howflow rates through narrow and intermediate-size nanotubes are related to bulk properties, and∗A version of this chapter has been published. L. Liu and G. N. Patey, J. Chem. Phys., 144, 184502(2016).149364.2. Water Flow Rates through a (6,6) CNTcan vary significantly for different water models, even when conduction occurs by the samemechanism.4.2 Water Flow Rates through a (6,6) CNTIn Figure 4.1 we plot typical cumulative counts of water molecules traveling through a(6,6) nanotube as functions of time. Because high pressures are applied in our simulations,molecules pass through the nanotube in a single direction from the high to low pressure reser-voir. Here we investigate water dynamics at four different temperatures: 260, 280, 300, and320 K. As for (8,8) and (9,9) nanotubes, the curves in Figure 4.1 are generally linear againsttime regardless of the temperature or the water model employed. As long as steady stateflows have been established, average flow rates can be easily obtained from linear fits of thecumulative counts curves. It is apparent that the flow rates are not identical for differentwater models, nor for the same model at different temperatures.The average flow rate is defined as the average numbers of water molecules traversingeither a (6,6) or a (9,9) nanotube during a period of 1 ns. The results of flow rates, along withestimated standard deviations are summarized in Table 4.1 for all three model (note that RCand RP models will be discussed in Chapter 4.4). The table includes flow rates measuredat four different temperatures as mentioned above, with the external pressure maintainedat 220 MPa. For both sizes of nanotube, the results demonstrate that the water transportrate can be strikingly influenced by the water model employed. If the size of the nanotubeand the temperature are fixed, TIP4P/2005 is always the slowest case, and TIP3P is alwaysthe fastest. The SPC/E model has intermediate flow rates, but in general, the SPC/E flowrates are closer to those of TIP4P/2005. Moreover, flow rates are positively correlated withtemperature for both (6,6) and (9,9) nanotubes. Interestingly, we find that ln(Rflow) has aclear linear dependence on T−1, where T is the absolute temperature, indicating that waterconduction is an activated process following an Arrhenius type equation374.2. Water Flow Rates through a (6,6) CNT0 2 4 6 8 10Time (ns)0100200300400Cumulative CountsTIP4P/2005SPC/ETIP3P(a)0 2 4 6 8 10Time (ns)0100200300400500Cumulative CountsTIP4P/2005SPC/ETIP3P(b)0 2 4 6 8 10Time (ns)0100200300400500600Cumulative CountsTIP4P/2005SPC/ETIP3P(c)0 2 4 6 8 10Time (ns)0200400600800Cumulative CountsTIP4P/2005SPC/ETIP3P(d)Figure 4.1: Cumulative counts curves showing the number of water molecules that passthrough a (9,9) nanotube as functions of time at different temperatures: 260 K (top left),280 K (top right), 300 K (bottom left), and 320 K (bottom right). Note that counts duringfirst 1 ns are not displayed because this is the equilibration period. The red, blue, and greencurves are for the TIP4P/2005, SPC/E, and TIP3P models, respectively.384.3. Explanation of the Model-dependent Water Transport Ratesln(Rflow) = − EaRT+ C , (4.1)where Ea is an activation energy, R is the gas constant, and C is effectively constant over thetemperature range considered. The calculated activation energies are also listed in Table 4.1.Additionally, for the TIP4P/2005 model we conducted simulations and calculated flow rateswith the external pressure fixed at 110 MPa, which is achieved by halving the external force fexin Equation (2.7). From Table 4.1, we note that the simulated flow rate is roughly proportionalto the external pressure, which is consistent with previous simulation results.140, 150, 151 Forthe (6,6) nanotube, the activation energies obtained at 220 and 110 MPa agree within theestimated standard deviations, implying that the barrier to flow is not strongly dependenton the driving pressure. For the (9,9) nanotube, the activation energy of TIP4P/2005 at 110MPa is significantly greater than that at 220 MPa, suggesting that the flow rate influencesthe barrier.4.3 Explanation of the Model-dependent Water TransportRatesIn Chapter 3 we discussed why water conduction through (8,8) and (9,9) CNTs is dependent onthe water model employed in simulations. Our argument was that the flow rate diversity stemsfrom different ring-like structures formed by the different water models inside the nanotubes.As a consequence, two distinctive conduction modes were proposed to explain the differentflow rates. However, in the (6,6) nanotube, the spatial restriction is such that each watermolecule can have at most two near neighbors, which results in the formation of a singlehydrogen-bonded chain. A snapshot of a typical molecular chain is shown in Figure 4.2. Thesingle hydrogen-bonded chain structure is a common feature of all three water models inthe (6,6) nanotube and gives rise to a so-called single-file conduction mode.47, 70, 95 This basicstructural feature and conduction mode are present at all temperatures considered. Therefore,unlike (8,8) and (9,9) cases, it is not possible to simply explain the different conduction rates394.3. Explanation of the Model-dependent Water Transport RatesPressure(MPa)Water modelTemperature Activationenergy(kJ mol−1)260 K 280 K 300 K 320 KFlow rates through a (6,6) CNT (ns−1)220TIP4P/2005 6 (1) 12 (1) 22 (3) 28 (2) 18.2 (2.0)SPC/E 12 (1) 21 (2) 34 (2) 42 (3) 14.8 (1.2)TIP3P 37 (3) 48 (2) 59 (2) 75 (3) 8.0 (1.0)RC 22 (1) 34 (1) 47 (2) 61 (3) 11.8 (0.7)RP 7 (1) 14 (1) 22 (2) 32 (2) 17.4 (1.7)110 TIP4P/2005 3 (1) 7 (2) 11 (1) 16 (3) 19.1 (4.2)Entry rates into a (6,6) CNT (ns−1)220TIP4P/2005 10 (1) 19 (1) 30 (3) 44 (1) 17.0 (1.2)SPC/E 18 (1) 30 (2) 45 (3) 59 (2) 13.8 (0.8)TIP3P 46 (2) 62 (3) 75 (2) 95 (2) 8.2 (0.5)RC 27 (2) 43 (1) 59 (2) 77 (4) 12.0 (1.0)RP 10 (1) 20 (1) 31 (1) 44 (2) 17.0 (1.2)110 TIP4P/2005 7 (1) 14 (2) 22 (2) 33 (1) 17.7 (2.8)Flow rates through a (9,9) CNT (ns−1)220TIP4P/2005 16 (4) 44 (4) 107 (6) 223 (7) 30.5 (2.7)SPC/E 31 (5) 67 (7) 134 (14) 292 (10) 25.6 (1.9)TIP3P 175 (10) 356 (6) 492 (15) 606 (11) 14.2 (0.7)110 TIP4P/2005 4 (2) 18 (4) 52 (5) 106 (3) 37.9 (5.5)Entry rates into a (9,9) CNT (ns−1)220TIP4P/2005 47 (6) 87 (4) 157 (7) 287 (7) 20.8 (1.4)SPC/E 66 (6) 108 (5) 187 (20) 359 (10) 19.3 (1.1)TIP3P 226 (11) 412 (6) 552 (10) 669 (8) 12.4 (0.5)110 TIP4P/2005 37 (8) 55 (5) 110 (9) 198 (3) 19.6 (2.4)Table 4.1: Summary of flow rates, entry rates, and their calculated activation energies fordifferent water models. The numbers in brackets are estimated standard deviations.404.3. Explanation of the Model-dependent Water Transport RatesFigure 4.2: A configurational snapshot of TIP4P/2005 water in the (6,6) nanotube. Carbonatoms are cyan, oxygen atoms are red, and hydrogen atoms are gray. The black dotted linesindicate hydrogen bonds. Note that part of carbon nanotube is not displayed for visual clarity.by invoking obvious structural differences within the (6,6) nanotube.To explain the mechanism of the model-dependent conduction and activation energy ob-served in the (6,6) case, we quantitatively examined properties of confined water that couldaccount for such differences. In particular, we calculated the average number of hydrogenbonds per water molecule, and the radial pressure exerted on the nanotube wall.Hydrogen bonds are defined using geometric criteria as described in Chapter 2.6.2, and fordifferent models and temperatures, the average numbers of hydrogen bonds per water molecule(hydrogen bond numbers) are tabulated in Table 4.2. The hydrogen bond number decreaseswith increasing temperature for all water models, implying that the chain structure becomesless rigid at higher temperatures. TIP3P has notably smaller hydrogen bond numbers thanthe other two models, which suggests that the degree of hydrogen bonding might be an aspectinfluencing the conduction rate. The TIP4P/2005 and SPC/E models have approximatelyidentical hydrogen bond numbers at all temperatures simulated, yet at 260 K, TIP4P/2005has half as many conduction counts as SPC/E (6 versus 12 ns−1), and the TIP4P/2005counts remain significantly smaller at higher temperatures. For this reason, the hydrogenbond number does not appear to determine the flow rate directly. From Table 4.2, we canalso conclude that the hydrogen bond number is in general independent of the magnitude ofthe external pressure.The pressure in the radial direction implicitly depends on water structural propertieswithin a CNT.152, 153 One would expect the radial pressure to influence the “friction” with414.3. Explanation of the Model-dependent Water Transport RatesPressure(MPa)Water modelTemperature260 K 280 K 300 K 320 KHydrogen bond number220TIP4P/2005 1.79 1.75 1.72 1.70SPC/E 1.78 1.76 1.73 1.71TIP3P 1.65 1.62 1.59 1.55RC 1.73 1.70 1.66 1.62RP 1.77 1.75 1.72 1.68110 TIP4P/2005 1.79 1.75 1.72 1.69Table 4.2: Average numbers of hydrogen bonds per water molecule (hydrogen bond numbers)within the (6,6) nanotube for different water models.the nanotube wall and hence possibly affect the flow rate. We estimated the radial pressure,Pradial, assuming that the equilibrium expression153Pradial =kBT 〈Nnt〉Vnt+12Vnt[〈∑i∑j>if(x)ij xij〉+〈∑i∑j>if(y)ij yij〉+〈∑if(i)nanotube(reff − ri)〉] (4.2)holds, where kB is the Boltzmann constant, T is the absolute temperature, 〈Nnt〉 is the averagenumber of water molecules within the nanotube, V is the volume of nanotube cavity, f(k)ij isthe k component of the force exerted by particle j on particle i, kij is the k component of thevector from j to i, f(i)nanotube is the force exerted by particle i on the nanotube wall, reff is theeffective radius of the nanotube, and ri is the distance of particle i from the symmetry axisof the nanotube.The calculated radial pressures are listed in Table 4.3. The pressure magnitudes arelarge and consistent with earlier equilibrium calculations.153 Given that for the TIP4P/2005model the radial pressure obtained does not have a strong dependence on the flow rate, theequilibrium assumption mentioned above appears reasonable. Also, the pressure exerted on424.4. Water Hydrodynamics and Bulk Fluid PropertiesPressure(MPa)Water modelTemperature260 K 280 K 300 K 320 KRadial pressure (MPa)220TIP4P/2005 205 (7) 208 (4) 216 (6) 227 (12)SPC/E 208 (5) 216 (3) 222 (4) 229 (7)TIP3P 179 (6) 190 (6) 199 (6) 204 (4)RC 202 (4) 205 (9) 217 (9) 224 (11)RP 213 (10) 216 (7) 226 (13) 237 (14)110 TIP4P/2005 208 (6) 218 (8) 228 (9) 233 (4)Table 4.3: Radial pressures exerted on the (6,6) nanotube by different water models. Thenumbers in brackets are estimated standard deviations.the nanotube wall increases with increasing temperature as expected. We note that the largestcontribution to the radial pressure (∼75%) comes from the LJ part of the water-nanotubeinteraction. The TIP3P radial pressures are obviously lower than those for TIP4P/2005 andSPC/E, possibly accounting for some of the faster flow rates. However, the radial pressures aresimilar between TIP4P/2005 and SPC/E, with the SPC/E values even slightly larger. Thusthe radial pressure can not account for why SPC/E has faster flow rates than TIP4P/2005.4.4 Water Hydrodynamics and Bulk Fluid PropertiesFailing to find a satisfactory explanation for the different flow rates considering only waterproperties within the (6,6) CNT, we asked whether the entry of water molecules into the nan-otube might be an important model-dependent factor. One will guess that the different watermodel employed determines the entry rate of water molecules if the radius of the nanotube isfixed. To investigate this idea we calculated water entry rates, Rentry, for different water mod-els. The method for calculating the entry rate is described in Chapter 2.6.1. The entry ratesare tabulated in Table 4.1. We see that for both (6,6) and (9,9) nanotubes, the entry ratesfor the three water models differ just as the flow rates do, with TIP3P being the fastest andTIP4P/2005 the slowest. Also, the entry rate becomes higher as the simulation temperature434.4. Water Hydrodynamics and Bulk Fluid Propertiesrises. Unlike the flow rate, which is generally proportional to the external pressure difference,the entry rate decreases to a lesser extent (only 30% change on average) when the externalpressure difference drops from 220 to 110 MPa.Moreover, we find that the temperature dependence of not only the flow rate but also of theentry rate follows an Arrhenius-like equation. Arrhenius plots for ln(Rflow) and ln(Rentry) areshown in Figures 4.3 and 4.4, respectively. We can draw several interesting conclusions fromthe estimated activation energies. For the (6,6) nanotube, the TIP4P/2005 model has largestactivation energies for both flow and entry rates while the TIP3P model has the smallestactivation energies. The same trends are also apparent in the case of the (9,9) nanotube.Most interestingly, in the (6,6) nanotube the activation energies associated with flow andentry agree within the standard deviations, regardless of the external pressure difference usedin simulations. However, in the (9,9) nanotube the activation energy similarity is not true,particularly for TIP4P/2005 and SPC/E, where the flow activation energies are notably higherthan the entry values.Based on the results in Table 4.1, we hypothesize that for the (6,6) case the different flowrates primarily come from different entry rates, and, therefore, the flow activation energy isessentially the corresponding entry activation energy. For the (9,9) case, this observation isnot correct. In Chapter 3, we showed that TIP4P/2005 and SPC/E are more likely to formpolygonal ring structures within the (9,9) CNT, in contrast with the TIP3P model. A usefulgraphical illustration is given in Figure 4.5, where we plot the rate ratios Rflow/Rentry asfunctions of temperature. We note that for the (6,6) nanotube, this ratio is approximatelyconstant within the temperature range investigated, which is consistent with nearly equal flowand entry activation energies. Whereas for the (9,9) nanotube, the ratios for all three modelsare dependent on temperature to a certain extend, as might be expected if the flow and entryactivation energies are different. Note that the profiles of TIP4P/2005 and SPC/E showstrong temperature dependence, while the profile of TIP3P only shows slight dependence.To identify the origin of the energy barrier to entry, we calculated the average potentialenergy of water (including both water-water and water-nanotube interactions) as a function444.4. Water Hydrodynamics and Bulk Fluid Properties0.003 0.0032 0.0034 0.0036 0.0038 0.004T-1 (K-1)2345ln(Rflow)TIP4P/2005SPC/ETIP3P(a)0.003 0.0032 0.0034 0.0036 0.0038 0.004T-1 (K-1)34567ln(Rflow)TIP4P/2005SPC/ETIP3P(b)0.003 0.0032 0.0034 0.0036 0.0038 0.004T-1 (K-1)2345ln(Rflow)TIP4P/2005RC modelRP model(c)Figure 4.3: The dependence of ln(Rflow) on T−1. Results for the (6,6) nanotube [panels (a)and (c)], and for the (9,9) nanotube [panel (b)] are shown. The water models are indicatedin the legends. Note that RC and RP refer to modified TIP4P/2005 models as described inthe text.454.4. Water Hydrodynamics and Bulk Fluid Properties0.003 0.0032 0.0034 0.0036 0.0038 0.004T-1 (K-1)2345ln(Rentry)TIP4P/2005SPC/ETIP3P(a)0.003 0.0032 0.0034 0.0036 0.0038 0.004T-1 (K-1)4567ln(Rentry)TIP4P/2005SPC/ETIP3P(b)0.003 0.0032 0.0034 0.0036 0.0038 0.004T-1 (K-1)2345ln(Rentry)TIP4P/2005RC modelRP model(c)Figure 4.4: The dependence of ln(Rentry) on T−1. Results for the (6,6) nanotube [panels (a)and (c)], and for the (9,9) nanotube [panel (b)] are shown. The water models are indicatedin the legends. Note that RC and RP refer to modified TIP4P/2005 models as described inthe text.464.4. Water Hydrodynamics and Bulk Fluid Properties260 280 300 320Temperature (K)00.20.40.60.81Rflow/RentryTIP4P/2005SPC/ETIP3PRC modelRP model(a)260 280 300 320Temperature (K)00.20.40.60.81Rflow/RentryTIP4P/2005SPC/ETIP3P(b)Figure 4.5: Temperature dependence of the ratio Rflow/Rentry for the (6,6) [panel (a)] and(9,9) [panel (b)] nanotubes. The water models are indicated in the legends. Note that RCand RP refer to modified TIP4P/2005 models as described in the text.474.4. Water Hydrodynamics and Bulk Fluid Propertiesof position along the z coordinate. Consider a cylinder centered on the symmetry axis ofthe (6,6) CNT, whose radius is equal to that of the nanotube, and which includes both thenanotube cavity and an extended region into the bulk. The cylinder is divided into bins each0.2 nm wide. Using equilibrium NV T simulations at 300 K, the average potential energy isdetermined for every bin. Corresponding potentials of mean force (PMF) for water along thez coordinate are also obtained as described in Chapter 2.6.4. Both thermodynamic profilesare obtained using five independent simulations each 100 ns long. The potential energy andthe PMF curves for all three water models are shown in Figures 4.6 and 4.7, respectively.We note that as a water molecule enters the nanotube, its interaction with other watermolecules decreases, but water-nanotube interactions can at least partially compensate sucheffect. From Figure 4.6 we note that the TIP3P model has almost no energy barrier to entry.The energy difference for TIP3P, defined as ∆U = UCNT − Ubulk, is ∼ −6.2 kJ mol−1. Thisobservation agrees well with the result (∼ −5.8 kJ mol−1) obtained by Waghe et al.53 forTIP3P water and a completely submerged (6,6) nanotube of similar length. The TIP4P/2005and SPC/E models have small energy barriers to entry (∼2.9 and ∼0.9 kJ mol−1). However,these values are not sufficient to account for the activation energies of Rentry (∼17.0 and ∼13.8kJ mol−1) even if we consider the substantial uncertainties at the nanotube orifice.The PMF profiles for all three water models shown in Figure 4.7 are qualitatively similar.However, the TIP4P/2005 and SPC/E models have considerable oscillation along the (6,6)nanotube symmetry axis, which we believe comes from increased hydrogen bonding withinthe nanotube (Table 4.2). The free energy barrier to entry for TIP3P we obtained (∼2.6 kJmol−1) is close to the value (∼2.1 kJ mol−1) reported by Corry.95 The free energy barriersare larger (∼3.0 kJ mol−1) for the TIP4P/2005 and SPC/E models, but again these energiesare too small to justify the observed differences in entry rate satisfactorily.Since both equilibrium potential energy and potential of mean force profiles fail to explainthe apparent barriers to entry, we focus our attention on another factor, the bulk diffusivityfor different water models. The method of determining self-diffusion coefficients has been de-scribed in Chapter 2.6.3. The diffusion coefficients, which agree well with literature values,132484.4. Water Hydrodynamics and Bulk Fluid Properties0 1 2 3 4Position (nm)-45-40-35Potential Energy (kJ mol-1 )Bulk region CNT regionTIP4P/2005SPC/ETIP3PFigure 4.6: Water potential energy profiles along the nanotube axis as described in the text.The water models are indicated in the legends. The error bars represent one standard devia-tion. The black vertical dotted line indicates the position of the (6,6) nanotube orifice.0 1 2 3 4Position (nm)12345PMF (kJ mol-1 )Bulk region CNT regionTIP4P/2005SPC/ETIP3PFigure 4.7: Water potential of mean force profiles along the nanotube axis as described in thetext. The water models are indicated in the legends. The black vertical dotted line indicatesthe position of the (6,6) nanotube orifice.494.4. Water Hydrodynamics and Bulk Fluid Properties0.003 0.0032 0.0034 0.0036 0.0038 0.004T-1 (K-1)-12-11-10-9ln(Dwater)TIP4P/2005 & RP modelSPC/ETIP3PRC modelFigure 4.8: The dependence of ln(Dwater) on T−1. Self-diffusion coefficients are in cm2 s−1.The water models are indicated in the legends.together with the associated activation energies are presented in Table 4.4. Arrhenius plotsfor the bulk diffusion coefficients are given in Figure 4.8. It is interesting to observe that theactivation energies of diffusion coincide in magnitude with those estimated from entry rates forall three water models, and display the same trend. TIP3P gives the lowest activation energyof diffusion (∼10.4 kJ mol−1), and TIP4P/2005 has the highest (∼19.6 kJ mol−1). Althoughthe activation energies obtained from entry rates and self-diffusion coefficients are not exactlyequal, they suggest that the entry rates are closely related to bulk fluid dynamics. Hence,even for single file conduction, different water models can feature different flow rates relatedto their different bulk mobilities. In the (6,6) nanotube case, the bulk diffusivity appears tobe much more important than any other factor in determining the flow rate. We supposethat we should take such factor into account in interpreting different flow rates obtained fromdifferent water models, and when attempting to compare simulation results with those of realwater.Several papers142, 154–156 have discussed entrance effect on nanotube flow rates. Mostinterestingly, Gravelle et al.156 demonstrated good agreement between MD simulation resultsand continuum hydrodynamics even for a nanotube undergoing single-file water conduction.504.4. Water Hydrodynamics and Bulk Fluid PropertiesWater modelTemperature Activation Energy(kJ mol−1)260 K 280 K 300 K 320 KSelf-diffusion coefficient (10−5 cm2s−1)TIP4P/2005and RP0.62 (0.03) 1.26 (0.02) 2.18 (0.10) 3.39 (0.04) 19.6 (0.6)SPC/E 0.88 (0.04) 1.63 (0.06) 2.62 (0.12) 3.67 (0.09) 16.5 (0.6)TIP3P 2.79 (0.07) 4.17 (0.07) 5.48 (0.25) 6.90 (0.06) 10.4 (0.3)RC 2.19 (0.07) 3.38 (0.09) 4.62 (0.21) 5.99 (0.29) 11.6 (0.6)Table 4.4: Water self-diffusion coefficients and their calculated activation energies for differentwater model. The numbers in brackets are estimated standard deviations.Therefore it is interesting to examine our data given continuum hydrodynamics.Sampson157 considered liquid flowing through a circular hole in an infinitely thin mem-brane and proposed the relationshipR =∆Pr3eff3η, (4.3)where R is either the fluid flow rate or entry rate, ∆P is the pressure difference, reff is theeffective radius of the hole, and η is the fluid shear viscosity. Note that in this case, theflow and entry rates are equivalent. Gravelle et al.156 even considered the hydrodynamicresistance of a cylindrical pore. They concluded that the rates are inversely proportional tothe shear viscosity. Also, using the criteria proposed by Gravelle et al.,156 the resistance ofthe nanotube in our specific case is negligible, therefore, Equation (4.3) remains valid. TheStokes-Einstein equation gives the connection between the diffusion coefficient, Dwater, andthe shear viscosity ηDwater =kBT6πηRStokes, (4.4)where RStokes is the Stokes radius (of water). From Equations (4.3) and (4.4), we easily obtainRTDwater=2πr3effRStokeskB∆P. (4.5)514.4. Water Hydrodynamics and Bulk Fluid PropertiesWater model σO−C (nm)εO−C(kJ mol−1)qH (e) qO (e) qM (e)TIP4P/2005 3.2793 0.5280 +0.5564 0 −1.1128RC 3.2793 0.5280 +0.5286 0 −1.0572RP 3.2793 0.5016 +0.5564 0 −1.1128Table 4.5: Force field parameters of the original TIP4P/2005, along with the modified RCand RP models.Thus, assuming that the RStokes is approximately the same for the different water models, theleft-hand side, RT/Dwater, should be constant regardless of the water model employed in MDsimulations.We plot the ratios RflowT/Dwater and RentryT/Dwater (in arbitrary units) against temper-ature in Figure 4.9. It is apparent that for both (6,6) and (9,9) CNTs the RentryT/Dwater,although slight differences amongst the models exist, remains nearly constant as the tem-perature varies for all models considered, which indicates that the entry rates are primarilydetermined by the hydrodynamics of the bulk liquid. For the (6,6) nanotube, RflowT/Dwaterbehaves much like RentryT/Dwater for all models, which is consistent with our conclusion thatthe entry rate largely decides the flow rate in this case. For the (9,9) nanotube, RflowT/Dwaterbehaves as RentryT/Dwater only for the TIP3P model. For the TIP4P/2005 and SPC/E mod-els, RflowT/Dwater differs with temperature and displays model dependence, which is againconsistent with our observation that water structures within the (9,9) nanotube have a stronginfluence on the flow rate.To further test the relationship between the flow rate and the bulk mobility, we conductedsimulations for two modified versions of the TIP4P/2005 model. In one modified model(denoted as RC), all charges in TIP4P/2005 water are reduced by 5%. In the other model(denoted as RP), the oxygen-carbon LJ energy parameter (εO−C) is reduced by 5%. Notethat in the RP model, only the water-nanotube interactions are altered, and the water-water interactions remain unchanged. The relevant force field parameters of the originalTIP4P/2005, RC, and RP models are given in Table 4.5.524.4. Water Hydrodynamics and Bulk Fluid Properties260 280 300 320Temperature (K)110100R flowT/DwaterTIP4P/2005SPC/ETIP3PRC modelRP model(a)260 280 300 320Temperature (K)110100RentryT/DwaterTIP4P/2005SPC/ETIP3PRC modelRP model(b)Figure 4.9: Temperature dependence of the ratio RflowT/Dwater [panel (a)] and RentryT/Dwater[panel (b)] in arbitrary units. Results for the (6,6) and (9,9) nanotubes are plotted as solid anddashed-dotted lines. The water models are indicated in the legends. Note that a logarithmicscale is used on the vertical axis.534.5. SummaryResults and analysis for the RC and RP models are included in Tables 4.1, 4.2, 4.3, 4.4,and in Figures 4.3, 4.4, 4.5, 4.8, 4.9. We observe that the RC model has faster flow and entryrates than those of the original TIP4P/2005 model, with correspondingly lower activationenergies. Besides, by reducing the partial charges on interaction sites and therefore weakeningthe overall water-water interactions, the bulk self-diffusion coefficients of RC are greater thanthat of TIP4P/2005. Moreover, the associated activation energy of RC is lower. Generallyspeaking, the results of the RC model are quite analogous to those of TIP3P, supportingour view that bulk transport essentially explains the different flow rates through the (6,6)CNT. Also, we see that the flow and entry rates of the RP model are almost identical tocorresponding rates of unmodified TIP4P/2005, even though the RP water interacts moreweakly with carbon atoms. This is also true of the activation energies. We remark that ourflow rate results agree with a previous study85 which mainly discussed the influence of thewater-nanotube interaction on the flow rate. Finally, although there are substantial differencesin the flow and entry rates between the two modified models, we see from Figure 4.9 that theresults again reasonably match the continuum hydrodynamic predictions.By investigating the dynamics of the two modified models, we demonstrate that the bulkfluid properties have a much greater impact on the water flow rate through the (6,6) nanotubethan details of the water-nanotube interaction.4.5 SummaryIn this chapter, we have shown that water conduction rates through a (6,6) CNT are de-pendent on the water model employed in MD simulations. This observation is consistentwith our earlier results for (8,8) and (9,9) nanotubes described in Chapter 3. However, un-like intermediate-size nanotubes in which confined water has special structures that influenceconduction rates, in the (6,6) nanotube all three models examined have similar single-stringconfigurations.We show that neither the different water structures within a (6,6) nanotube nor different544.5. Summarywater-carbon interactions can account for the model dependence of the flow rate. We discoverthat both the water flow rate through and the entry rate into the (6,6) nanotube are activatedprocesses. Moreover, the temperature dependences of both rates are well described by anArrhenius-like equation. Based on estimated activation energies, we reason that the flowrates are closely related to the entry rates, which in turn are strongly influenced by the bulkmobilities of the different water models. Continuum hydrodynamics calculations, as well asMD simulation results employing modified water models, support this conclusion.Our results unequivocally demonstrate that the water flow rate through even a (6,6)nanotube can be strikingly model dependent. We trace this somewhat unexpected observationto the different diffusion abilities of bulk water. Given the present observations, it would bedesirable to select a water model such as TIP4P/2005, which has a self-diffusion coefficientclose to that of real water, when attempting to simulate water transport through CNTs orother nanoscopic channels.55Chapter 5How Different Water Models AffectSimulated Ion Transport Ratesthrough a (9,9) CNT∗5.1 OverviewIn previous chapters, we showed that different water models can have strikingly different flowrates through CNTs. To be more specific, for intermediate-size nanotubes such as the (9,9)nanotube, not only the bulk mobility but also the structure of water confined within thenanotube determines the flow rate. In this chapter, we consider NaCl solutions with con-centrations varying from ∼0.25 to ∼2.8 mol L−1 and investigate water and ion conductionthrough a finite-length (9,9) nanotube. We observe that water and ion flow rates still differsharply for the TIP4P/2005 and TIP3P models. We note that the flow rate dependence ontemperature fits an Arrhenius-type equation. Also, for the TIP3P model, simulated watertransport rates in solution agree well with expectations based on continuum hydrodynamics.Most importantly, we confirm that both factors, bulk fluid diffusivity and confined molecularstructure, again can account for ion transport differences between models. Our results demon-strate that ion conduction through nanotubes is sensitive to factors other than ion itself, inparticular, the structure of the confined water molecules.∗A version of this chapter has been published. L. Liu and G. N. Patey, J. Chem. Phys., 146, 074502(2017).158565.2. Water and Ion Flow Rates through a (9,9) CNT5.2 Water and Ion Flow Rates through a (9,9) CNTHere we focus on only two water models: TIP4P/2005 and TIP3P. In Chapters 3 and 4 weshowed that these two models have the most significant differences. Although the TIP4P/2005model has been rarely studied in simulations of water transport, it is probably the most real-istic rigid water model available.132 The TIP3P model, which has less accurate bulk transportproperties, has been more widely used in water transport studies.88, 90, 95, 107, 111, 151, 159In Figure 5.1 we plot typical cumulative counts for water molecules and ions in NaClsolution passing through a (9,9) nanotube as functions of time. The results shown are for a1 mol L−1 solution at 320 K, and we note that the counts curves are generally linear againsttime for both water models, and at all solution concentrations and temperatures. Note thatfor ions, especially the TIP4P/2005 solution, the flow is rather slow resulting in jagged curves,nevertheless, an increasing trend is clear with essentially a constant rate. The average flowrates, Rflow, for both water and ions can be estimated by linearly fitting the cumulativecounts curves after a steady state is completely established. These rates, as well as standarddeviations, are tabulated in Table 5.1.We first considered water flow rates. We see immediately from Table 5.1 that the water flowrate significantly depends on the water model employed, with TIP3P having much faster flowrates than TIP4P/2005 at the same temperature and salt concentration. This is especiallytrue at the lower temperatures. With increasing temperature, faster water flow rates areobserved for both models. This observation is not surprising given our earlier discussionsin Chapters 3 and 4, showing that both the fluid viscosity and water structure inside thenanotube can play important roles in determining water conduction. Moreover, the waterflow rate generally decreases with ion concentration by much more than the small decrease inthe number of water molecules in the MD simulations. We believe that the decrease can be atleast partially explained by the decreasing diffusion coefficients of bulk water with increasingsalt concentration. For example, at 300 K the diffusion coefficients of TIP4P/2005 water are2.18×10−5, 2.05×10−5, 1.73×10−5, and 1.10×10−5 cm2 s−1, for pure water and NaCl solutions575.2. Water and Ion Flow Rates through a (9,9) CNT0 2 4 6 8 10Time (ns)010002000300040005000Cumulative CountsTIP4P/2005TIP3P(a)0 2 4 6 8 10Time (ns)010203040Cumulative CountsTIP4P/2005, Na+TIP4P/2005, Cl-TIP3P, Na+TIP3P, Cl-(b)Figure 5.1: Cumulative counts curves showing the number of water molecules [panel (a)] andions [panel (b)] that pass through a (9,9) nanotube as functions of time. Results for water areshown as in solid curves, and sodium and chloride ions are represented by dashed and dottedcurves, respectively. Note that counts during first 1 ns are not displayed because this is theequilibration period. The red and green curves are for the TIP4P/2005 and TIP3P models,respectively.585.2. Water and Ion Flow Rates through a (9,9) CNTWatermodelParticleTemperature Activationenergy(kJ mol−1)260 K 280 K 300 K 320 KFlow rates (ns−1), 0.25 mol L−1 NaCl solutionTIP4P/2005H2O 10.1 (2.8) 37.3 (4.4) 102.1 (2.9) 206.3 (8.4) 35.0 (3.1)Na+ < 0.1 < 0.1 0.1 (0.2) 0.3 (0.1)n/aCl− < 0.1 < 0.1 0.1 (0.1) 0.2 (0.2)TIP3PH2O 174.6 (6.6) 331.0 (9.2) 452.9 (11.9) 570.9 (9.0) 13.5 (0.5)Na+ 0.4 (0.2) 0.8 (0.4) 1.1 (0.3) 1.4 (0.4) 14.3 (6.4)Cl− 0.2 (0.1) 0.4 (0.2) 0.5 (0.3) 0.6 (0.2) 12.4 (6.9)Flow rates (ns−1), 1 mol L−1 NaCl solutionTIP4P/2005H2O 8.6 (4.2) 34.3 (3.0) 89.8 (5.9) 179.8 (5.6) 35.1 (5.3)Na+ < 0.1 < 0.1 0.4 (0.2) 1.3 (0.4)n/aCl− < 0.1 < 0.1 0.2 (0.2) 0.8 (0.1)TIP3PH2O 175.6 (6.8) 283.5 (17.9) 389.7 (7.6) 488.0 (9.4) 11.8 (0.5)Na+ 1.3 (0.2) 2.3 (0.3) 3.1 (0.4) 4.1 (0.6) 13.1 (2.3)Cl− 1.0 (0.4) 1.7 (0.5) 2.2 (0.5) 2.9 (0.2) 12.0 (4.6)Flow rates (ns−1), 2.8 mol L−1 NaCl solutionTIP4P/2005H2O 7.4 (1.5) 27.2 (7.5) 81.3 (4.1) 145.2 (7.3) 34.9 (2.4)Na+ < 0.1 0.2 (0.2) 1.3 (0.2) 2.8 (0.5)n/aCl− < 0.1 0.2 (0.2) 1.0 (0.3) 2.2 (0.3)TIP3PH2O 108.5 (15.9) 196.4 (6.4) 297.5 (10.8) 372.8 (10.6) 14.4 (1.6)Na+ 2.0 (0.6) 3.9 (0.6) 5.5 (0.5) 7.0 (0.7) 14.4 (3.4)Cl− 1.8 (0.7) 3.0 (0.7) 4.8 (0.4) 5.9 (0.8) 14.0 (4.5)Table 5.1: Summary of flow rates for water molecules and ions as well as the calculatedactivation energies for different water models. The numbers in brackets are estimated standarddeviations.595.2. Water and Ion Flow Rates through a (9,9) CNT0 0.5 1 1.5 2 2.5 3Concentration (mol L-1)02468R flow,ion (ns-1 )TIP4P/2005, Na+TIP4P/2005, Cl-TIP3P, Na+TIP3P, Cl-Figure 5.2: Ion flow rates as functions of the NaCl concentration. Results are obtained fromsimulations conducted at 320 K. The water models and specific particles are indicated in thelegend.at 0.25, 1, and 2.8 mol L−1, respectively. The corresponding diffusion coefficients for TIP3Pare 5.48×10−5, 5.35×10−5, 4.64×10−5, and 3.23×10−5 cm2 s−1, displaying the same trend.While water flow rates remain measurable at low temperatures, the ion flow rates aresometimes too slow to be determined in 10 ns simulations, especially for the TIP4P/2005model. For example, at 260 K it required on average 20 ns in order to observe a single ionconduction event for the TIP4P/2005 model. In general, ion flow rates for TIP3P solutionsare much faster than those for TIP4P/2005. We also note that Na+ has faster transportrates than Cl− for both water models. The ion flow rate is also sensitive to temperature,increasing as the temperature increases. Additionally, the ion flow rate increases with the saltconcentration of the feed reservoir, as illustrated in Figure 5.2.We would expect the ion flow rate to be strongly influenced by the water flow rate and theion concentration. Therefore, in order to control for these influences and isolate other factorsaffecting ion conduction, it is useful to introduce an ion transport efficiency parameter, ξion,defined as605.3. Why Water and Ion Conduction Is Dependent on the Water Model Employedξion =Rflow,ion/ρionRflow,water/ρwater, (5.1)where ρwater and ρion are particle number densities for water molecules and either ion, re-spectively. Note that if the ion and water flow rates normalized by the particle densitieswere equal, then ion transport efficiency parameters would be one. Parameters less than oneindicate that ion transport is less efficient than water transport. It is worth mentioning thata somewhat similar quantity, called “the percentage salt rejection”, was defined by Corryand co-workers88, 95, 159 to quantify the resistance to ion conduction through a membraneconstructed of nanotubes.Results for the ion transport efficiency parameter from MD simulations are summarizedin Table 5.2. For both water models, the parameters are always significantly less than one,suggesting that ions experience more hindrance than water during transport. Also, the sodiumion parameters are larger than or equal to those of chloride ion, as is obvious from the relativeion flow rate values given in Table 5.1. However, apart from these common features, theparameters for TIP4P/2005 and TIP3P solutions are strikingly different. At 260 and 280 K,the parameters of both sodium and chloride ions for the TIP4P/2005 model are too smallto accurately measure even in our longest simulation runs, whereas, for the TIP3P modelthe parameters are significant and readily accessible at these temperatures. Also, the iontransport efficiency parameter is approximately constant for TIP3P but strongly influencedby temperature for TIP4P/2005.5.3 Why Water and Ion Conduction Is Dependent on theWater Model EmployedIn earlier chapters, we demonstrated that for the (9,9) CNT, pure water transport is dependenton two factors: the bulk fluid shear viscosity and the structure of confined water. Based onthese observations, a logical speculation is that the model dependence of the ion flow rateis of similar origin. We argue that there are two potential explanations. It is possible that615.3. Why Water and Ion Conduction Is Dependent on the Water Model EmployedWater model ParticleTemperature260 K 280 K 300 K 320 KIon transport efficiency parameter, 0.25 mol L−1 NaCl solutionTIP4P/2005Na+n/a n/a0.22 (0.42) 0.32 (0.09)Cl− 0.22 (0.21) 0.21 (0.20)TIP3PNa+ 0.50 (0.23) 0.53 (0.25) 0.53 (0.13) 0.54 (0.15)Cl− 0.25 (0.12) 0.26 (0.19) 0.24 (0.14) 0.23 (0.07)Ion transport efficiency parameter, 1 mol L−1 NaCl solutionTIP4P/2005Na+n/a n/a0.24 (0.10) 0.39 (0.11)Cl− 0.18 (0.11) 0.24 (0.02)TIP3PNa+ 0.40 (0.05) 0.43 (0.03) 0.42 (0.05) 0.45 (0.03)Cl− 0.30 (0.08) 0.32 (0.07) 0.30 (0.06) 0.32 (0.02)Ion transport efficiency parameter, 2.8 mol L−1 NaCl solutionTIP4P/2005Na+n/a0.13 (0.10) 0.29 (0.03) 0.35 (0.04)Cl− 0.13 (0.16) 0.22 (0.06) 0.27 (0.02)TIP3PNa+ 0.33 (0.05) 0.36 (0.04) 0.33 (0.02) 0.34 (0.02)Cl− 0.30 (0.07) 0.28 (0.06) 0.29 (0.01) 0.29 (0.03)Table 5.2: Ion transport efficiency parameters ξion for different water models. The numbersin brackets are estimated standard deviations.625.3. Why Water and Ion Conduction Is Dependent on the Water Model Employedthe ion flow rate is closely connected with the mobility of the bulk solution, and the ring-like structures occurring with TIP4P/2005 impede not only water but also ion conduction.To further explore the mechanism of ion conduction and its dependence on the water modelemployed, in the following analysis we focus on the 1 mol L−1 system.We find that flow rate results of water again follow the Arrhenius relationship (Equation(4.1)), as plotted in Figure 5.3. The estimated activation energies are included in Table 5.1.As with pure water, TIP4P/2005 has a much higher activation energy than TIP3P (∼35.1versus ∼11.8 kJ mol−1) in salt solutions. For both TIP4P/2005 and TIP3P, the activationenergies do not strongly depend on the solution concentration. Note that the activationenergy for TIP4P/2005 obtained from 1 mol L−1 simulations is higher than that of purewater simulations (∼35.1 versus ∼30.5 kJ mol−1). We believe that both activation energiesagree with each other if we allow for the large standard deviations of activation energies (5.3and 2.7 kJ mol−1). The ion flow rates for TIP3P are well fit by the Arrhenius relationship,as plotted in Figure 5.4. But for the TIP4P/2005 model the very low flow rates of ions andassociated large uncertainties at 260 and 280 K do not allow meaningful fits to the Arrheniusequation.Entry rate results, Rentry, together with activation energies are presented in Table 5.3.We observe that for TIP4P/2005, the estimated activation energy of the water entry rateobtained from the 1 mol L−1 solution agrees with that from the pure water system within theerror estimates. Whereas for TIP3P, such activation energies are close with a difference of∼2.0 kJ mol−1. For TIP4P/2005 ion entry rates are not too small to be accurately measuredin our simulations. It is interesting to compare ion entry rates with ion flow rates. For TIP3Pthe entry rates for chloride ions are faster or, taking account of error estimates, at least equalto those of sodium ions, whereas the chloride flow rates are lower than those of sodium. Forthe TIP4P/2005 model, the entry rates for both ions are indistinguishable within the errorestimates, but, at least at 320 K, the sodium ion conduction rate is faster than the chlorideion. It suggests that sodium ions experience less resistance than chloride ions, which is likelydue to the smaller “size” of the sodium ion.160635.3. Why Water and Ion Conduction Is Dependent on the Water Model Employed0.003 0.0032 0.0034 0.0036 0.0038 0.004T-1 (K-1)234567ln(Rflow,water)TIP4P/2005TIP3PFigure 5.3: The dependence of ln(Rflow,water) on T−1. The results are obtained from simula-tions of 1 mol L−1 NaCl solution; results at other concentrations show similar behavior. Thewater models are indicated in the legends.0.003 0.0032 0.0034 0.0036 0.0038 0.004T-1 (K-1)012ln(Rflow,ion)TIP3P, Na+TIP3P, Cl-Figure 5.4: The dependence of ln(Rflow,ion) on T−1. The results are obtained from simulationsof 1 mol L−1 NaCl solution; results at other concentrations show similar behavior. The watermodel and specific particles are indicated in the legends.645.3. Why Water and Ion Conduction Is Dependent on the Water Model EmployedDiffusion coefficients of all three particles are summarized results in Table 5.4. The re-sults obtained for sodium and chloride ions with the TIP3P model at 300 K are close tothose (2.28×10−5 cm2 s−1 for sodium and 2.90×10−5 for chloride) reported by Joung andCheatham.137 Joung and Cheatham137 also pointed out the large discrepancies in the iondiffusion coefficients obtained with different water models. The activation energies includedin Table 5.4 are calculated from fits to the Arrhenius equation. Activation energies of iondiffusion coefficients are relatively close to those of ion entry rates, demonstrating that againthe entry rates are strongly influenced by the bulk viscosity.Plots of the rate ratios Rflow/Rentry for water and ions are given in Figure 5.5. For water,the pattern is the same as in our pure water simulations. The ratio for TIP3P is ∼0.8 andshows little dependence on temperature, suggesting that the flow and entry rates are closelyconnected at all temperatures investigated. In contrast, the ratio for TIP4P/2005 is ∼0.2 at260 K and increases to ∼0.7 at 320 K. For TIP4P/2005, we believe that the entry and flowrates do not strongly correlate at lower temperatures, where the water structure within thenanotube slows down the conduction, but the importance of the entry rate increases at highertemperatures where confined water is less structured. To confirm this reasoning, we againanalyze the water structure inside the CNT as described in Chapter 2.6.2. In particular, itis useful to determine the average percentage of water molecules within the nanotube thatare not part of any ring-like configuration, or ring-free molecules. The average number ofhydrogen bonds per water molecule (hydrogen bond number) is useful as well given thatring-free molecules often associate with a smaller hydrogen bond number. The results givenin Table 5.5 are obtained from MD simulations of 1 mol L−1 NaCl solution. Note that weused only frames where no ions were present within the nanotube in calculating the averages.From Table 5.5 we see that there are far fewer ring-free molecules for TIP4P/2005 than forTIP3P, indeed for TIP4P/2005 ring free molecules are very rare at the lower temperatures.For both water models, the percentage of ring-free molecules increases with temperature, butthe growth is much more dramatic for TIP4P/2005 (2.5% at 260 K to 41.9% at 320 K) thanfor TIP3P (32.5% to 65.6%, correspondingly). Consistently, hydrogen bond numbers for the655.3. Why Water and Ion Conduction Is Dependent on the Water Model EmployedWatermodelParticleTemperature Activationenergy (kJmol−1)260 K 280 K 300 K 320 KEntry rates (ns−1), 1 mol L−1 NaCl solutionTIP4P/2005H2O 44.7 (2.7) 84.6 (4.2) 153.5 (6.4) 257.3 (5.4) 20.2 (0.7)Na+ < 0.1 0.3 (0.3) 1.0 (0.2) 2.3 (0.5)n/aCl− < 0.1 0.2 (0.1) 0.8 (0.4) 2.4 (0.4)TIP3PH2O 235.9 (3.8) 355.1 (13.9) 470.6 (5.3) 579.7 (4.9) 10.4 (0.2)Na+ 2.0 (0.2) 3.1 (0.4) 3.9 (0.3) 4.9 (0.6) 10.2 (1.7)Cl− 2.2 (0.3) 3.5 (0.4) 4.1 (0.3) 5.2 (0.3) 9.6 (1.6)Table 5.3: Summary of entry rates for water molecules and ions in addition to the calculatedactivation energies. The numbers in brackets are estimated standard deviations.WatermodelParticleTemperature Activationenergy(kJ mol−1)260 K 280 K 300 K 320 KDiffusion coefficient (10−5 cm2·s−1), 1 mol L−1 NaCl solutionTIP4P/2005H2O 0.43 (0.02) 0.95 (0.04) 1.73 (0.07) 2.67 (0.08) 21.1 (0.6)Na+ 0.11 (0.01) 0.38 (0.02) 0.69 (0.07) 1.16 (0.17) 26.8 (1.8)Cl− 0.27 (0.08) 0.53 (0.05) 0.88 (0.07) 1.36 (0.46) 18.6 (4.7)TIP3PH2O 2.29 (0.06) 3.33 (0.16) 4.64 (0.16) 5.99 (0.10) 11.1 (0.4)Na+ 0.80 (0.27) 1.51 (0.28) 2.18 (0.57) 2.73 (0.29) 14.2 (4.0)Cl− 1.38 (0.54) 1.76 (0.26) 2.74 (0.59) 3.12 (0.49) 10.0 (4.6)Table 5.4: Summary of diffusion coefficients for water molecules and ions together with thecalculated activation energies. The numbers in brackets are estimated standard deviations.665.3. Why Water and Ion Conduction Is Dependent on the Water Model Employed260 280 300 320Temperature (K)00.20.40.60.81Rflow/RentryTIP4P/2005TIP3P(a)260 280 300 320Temperature (K)00.20.40.60.81Rflow/RentryTIP4P/2005, Na+TIP4P/2005, Cl-TIP3P, Na+TIP3P, Cl-(b)Figure 5.5: Temperature dependence of the ratio Rflow/Rentry for water [panel (a)] and ions[panel (b)]. The water models and specific particles are indicated in the legends.675.3. Why Water and Ion Conduction Is Dependent on the Water Model EmployedWater modelTemperature260 K 280 K 300 K 320 KRing-free molecules (%)TIP4P/2005 2.5 5.3 14.8 41.9TIP3P 32.5 57.4 61.8 65.6Hydrogen bond numberTIP4P/2005 3.39 3.33 3.22 3.01TIP3P 3.26 2.88 2.75 2.66Table 5.5: Selected structural properties of water molecules confined within the nanotube.TIP4P/2005 model are larger at all four temperatures.Considering Rflow/Rentry for ions shown in Figure 5.5 we see a similar picture. For TIP3Pthe ratios for sodium ions vary from ∼0.65 to ∼0.84 as the temperature rises from 260 to320 K, while the ratios of chloride ion over the same temperature range vary from ∼0.45 to∼0.56, which are remarkably smaller. We attribute this observation to the relative ion size.160The larger size of the chloride somewhat impedes its transport through the nanotube. ForTIP4P/2005 ratios of both ions are too small to be determined at 260 and 280 K, but at 300and 320 K, ratios for both ions are noticeably smaller than the TIP3P results, consistent withwater structure acting to reduce ion conduction.By inspecting configurational snapshots that show the ion environment within the CNTwe can gain additional insight on this issue. Typical snapshots of 1 mol L−1 systems at300 K are displayed in Figure 5.6. We notice that an ion is usually surrounded by ring-freemolecules within the nanotube. The snapshots also support our previous results (Table 5.5)which show that TIP4P/2005 has fewer ring-free molecules than TIP3P. Also, we see that forTIP4P/2005, the fraction of ring-free molecules and the ion transport efficiency parametersare both strongly temperature dependent, however, for TIP3P the temperature dependenceis much weaker. These results suggest that an ion can more easily enter the nanotube if themolecules near the orifice are ring-free, or in other words if any ring-like structures are brokendown before or during the entry process. Perhaps this is also true for ion passage through the685.3. Why Water and Ion Conduction Is Dependent on the Water Model Employednanotube.Equilibrium ion potential of mean force (PMF) curves for 1 mol L−1 systems at 300 Kare plotted in Figure 5.7. The PMF results provide some support for the argument givenabove. For the TIP3P model, the free energy barriers for sodium and chloride ions are ∼9.9and ∼10.6 kJ mol−1, which are in accord with earlier results obtained by Beu.97 We notethat the higher free energy barrier for chloride agrees well with its slightly smaller flow ratecompared with sodium, although the chloride ion does have a somewhat faster entry rate.For the TIP4P/2005 model, the equilibrium free energy barriers are obviously larger, ∼15 kJmol−1 for sodium and ∼20 kJ mol−1 for chloride ion, consistent with the much slower entryand flow rates observed for this water model. Moreover, the PMF curves for the TIP4P/2005solution increase as the ion moves towards the center of the nanotube, possibly suggestingadditional friction during the conduction process.We have already discussed entrance effects on water conduction through CNTs in Chapter4.4. It is natural to question whether such effects influence ion conduction as well. It has beenshown that continuum hydrodynamics can, at least in some cases, give a good description ofwater entry into (and sometimes flow through) nanotubes. For the (9,9) nanotube, we showedthat the ratio RT/Dwater of the TIP3P model is approximately constant for both entry andflow rates. For the TIP4P/2005 model, the relationship is roughly correct for the entry ratebut does not hold for the flow rate. We attributed this failure to the more complex molecularstructure formed within the (9,9) nanotube, which creates an additional resistance to flow.Here we considered RT/Dwater for both water and ions in the 1 mol L−1 solution. Providedthat RStokes is not strongly temperature dependent (which is reasonably anticipated), wewould expect the ratio to be approximately constant regardless of the simulated temperature.Plots of Rflow,waterT/Dwater and Rentry,waterT/Dwater are given in Figure 5.8. We notethat Rentry,waterT/Dwater shows little variation with temperature for both water models.Rflow,waterT/Dwater for TIP3P is also generally constant, but not for TIP4P/2005 where theratio increases with temperature. We see that the behavior pattern is much as we previouslyobserved for pure water for 1 mol L−1 solutions.695.3. Why Water and Ion Conduction Is Dependent on the Water Model EmployedFigure 5.6: Configurational snapshots of the ion environment with TIP4P/2005 (top) andTIP3P (bottom) water inside the nanotube at 300 K. Carbon atoms are cyan, the sodiumion is blue, the stacked ring structures are repeating green, orange and yellow, and ring-freewater molecules are magenta. Note that part of the carbon nanotube is not displayed forvisual clarity.0 1 2 3 4Position (nm)-100102030PMF (kJ mol-1 )Bulk region CNT regionTIP4P/2005 Na+TIP4P/2005, Cl-TIP3P, Na+TIP3P, Cl-Figure 5.7: Sodium and chloride ion potential of mean force profiles along the nanotube axisas described in the text. The water models and specific particles are indicated in the legends.The black vertical dotted line indicates the position of the (9,9) nanotube orifice.705.3. Why Water and Ion Conduction Is Dependent on the Water Model Employed260 280 300 320Temperature (K)102103104105R flow,waterT/DwaterTIP4P/2005, waterTIP3P, water0200400Rflow,ionT/DwaterTIP4P/2005, Na+TIP4P/2005, Cl-TIP3P, Na+TIP3P, Cl-(a)260 280 300 320Temperature (K)103104105R entry,waterT/DwaterTIP4P/2005, waterTIP3P, water0200400R entry,ionT/DwaterTIP4P/2005, Na+TIP4P/2005, Cl-TIP3P, Na+TIP3P, Cl-(b)Figure 5.8: Temperature dependence of the ratio RflowT/Dwater [panel (a)] and RentryT/Dwater[panel (b)] in arbitrary units. The water models and specific particles are indicated in thelegends. The scale on the left-hand y axis is for water, and the scale on the right-hand y axisis for ions.715.4. SummaryPlots of Rflow,ionT/Dwater and Rentry,ionT/Dwater are also included in Figure 5.8. We notethat for the TIP3P model, although Rflow,ionT/Dwater does appear to show a small increasewith temperature, Rentry,ionT/Dwater is again almost constant from 260 to 320 K. The ratioof sodium ions is greater than that of chloride ions. It is reasonable given that a hydratedsodium ion is larger than a hydrated chloride ion.160 For the TIP4P/2005 model, however,neither ratio remains constant as the temperature increases, indicating that for this modelthe hydrodynamic description is not adequate even for the entry rate. The ratio values givelarge uncertainties (not shown in Figure 5.8). One possible explanation of the fact thatRentry,ionT/Dwater has a strong temperature dependence is that the number of ring-free watermolecules near the orifice of the nanotube increases with increasing temperature, thus reducinga possible barrier to entry. We would not expect continuum hydrodynamics to capture suchstructural effects. It is worth noting that if the ratio is calculated using Dion rather thanDwater, similar plots are obtained.5.4 SummaryWe have examined pressure-driven conduction of water and ions through a (9,9) CNT. Resultsare presented for two common water models, TIP3P and TIP4P/2005, with NaCl concentra-tions ranging from 0.25 to 2.8 mol L−1, and temperatures from 260 to 320 K. Not surprisingly,our results for water conduction closely parallel our previous findings for pure water transport.For the TIP4P/2005 model, the entry rates are well described by continuum hydrodynam-ics, but the much slower flow rates are significantly impacted by the water structure withinthe nanotube, especially at lower temperatures. While for the TIP3P model, both waterflow and entry rates are highly connected, which can be largely understood by continuumhydrodynamics.Ion conduction exhibits similar model dependence. We find that the ion transport ef-ficiency parameters are much smaller for TIP4P/2005 than for TIP3P. In other words, ionconduction experiences stronger resistance with TIP4P/2005 solutions than with TIP3P. It725.4. Summaryis particularly true at lower temperatures where the flow rates are too low to be accuratelymeasured in our simulations. We provide evidence that the extra resistance results fromthe enhanced water structure that occurs with the TIP4P/2005 model within the nanotube.We note that the ion entry rates for the TIP3P model are approximately consistent withcontinuum hydrodynamics, but that is not true for TIP4P/2005 solutions.Our results have shown that water structure acts to impede the ion transport, while stillallowing water to flow at a significant rate. It is important because, as we demonstrate, itcan result in strong model dependence, and this should be kept in mind whenever molecularsimulations are being used to investigate ion flow through nanoscopic channels. Addition-ally, assuming that TIP4P/2005 is a reasonably good water model, our results indicate thatreducing the temperature will increase water structure, and strongly reduce ion transportefficiency.73Chapter 6Conclusions and Perspective6.1 Simulation Results of Water and Ion Transport throughCNTsWith the help of MD simulations, this dissertation investigates water and ion conductionthrough CNTs under external pressures. We are particularly interested in the influence ofwater models employed in simulations. Surprisingly, we have shown that the flow rates ofboth water and ions are remarkably sensitive to the water model employed through bothnarrow and intermediate-size nanotubes. To seek explanations, we conducted an analysis,investigating the bulk fluid dynamics, the confined water structure, and other potentiallyrelevant factors.In Chapter 3, we focus on water conduction across intermediate-size CNTs, specifically(8,8) and (9,9) nanotubes, which allow water molecules to form certain ordered arrangements,for example, stacked polygonal configurations, within their cavities. We find that differentwater models can have distinctive structural features, which can influence the mode of waterconduction. Generally speaking, a water model which forms ring-like configurations (rep-resented by TIP4P/2005) tends to adopt a slower cluster-by-cluster conduction mode withmany ring clusters moving as single units through the nanotube. In contrast, a water modelwhich forms less ordered structure (represented by TIP3P) is likely to choose a faster diffusiveconduction mode with molecules moving as separate particles, rather than as parts of largerclusters.Within narrow CNTs, such as a (6,6) nanotube, water forms only single-strand, hydrogen-746.1. Simulation Results of Water and Ion Transport through CNTsbonded chains for all three water models considered. Because characteristic structural fea-tures, which are observed in intermediate-size nanotubes, do not occur in the (6,6) case, ourprevious arguments explaining the different water flow rates do not hold. Thus, in Chapter 4we look for an explanation which can account for the model dependence observed for a narrownanotube. We discover that the water flow rate strongly correlates with the water entry rate.Both the water entry rate into, and flow rate through, CNTs are activated processes. AnArrhenius-type equation well describes the temperature dependences of both rates. Also, theentry rate appears closely related to the self-diffusion coefficient (or equivalently the shearviscosity) of the particular water model employed.For the (6,6) CNT, both the flow and entry rates for all water models examined areapproximately inversely proportional to the shear viscosity of the bulk liquid, as predictedby continuum hydrodynamics. Therefore, the differences in bulk mobility are reflected in thedifferent flow rates, and this accounts for the strong model dependence. For the (9,9) CNT,these observations also apply to a water model (represented by TIP3P) which has less orderedwater structure within the nanotube. However, although the entry rate again agrees withthe continuum hydrodynamics predictions for a water model (represented by TIP4P/2005)which has more ordered water structure within the nanotube, this does not apply to the flowrate. In all likelihood, the flow rate is strongly influenced by the water structure within theintermediate-size nanotubes.In Chapter 5 we investigate water and ion transport through a (9,9) CNT using NaClsolutions of different concentrations. We observe that the flow rates of water, sodium, andchloride ions through the nanotube are not only strongly model dependent but temperaturedependent as well. The water conduction rates can be explained as in the pure water simula-tions. We examine ion conduction in view of the two factors noted above. We show that theion mobilities in solution can have some influence on the ion conduction rate. In addition,we conclude that highly hydrogen-bonded water structure within the nanotube can greatlyincrease the resistance to ion conduction. Our results demonstrate that increasing the waterstructure within the CNT by decreasing the temperature strongly inhibits ion conduction,756.1. Simulation Results of Water and Ion Transport through CNTsWater modelTemperature Activationenergy(kJ mol−1)278 K 298 K 318 K 373 KSelf-diffusion coefficient (10−5 cm2s−1)Experimental 1.31 2.30 3.57 n/a 18.4TIP4P/2005 1.27 2.06 3.07 n/a 16.2SPC/E 1.54 2.54 3.57 n/a 15.4TIP3P 3.72 5.49 6.31 n/a 9.7Shear viscosity, 1 bar (10−4 Pa s)Experimental n/a 8.96 n/a 2.84 n/aTIP4P/2005 n/a 8.55 n/a 2.89 n/aSPC/E n/a 7.29 n/a 2.69 n/aTIP3P n/a 3.21 n/a 1.65 n/aTable 6.1: Experimental data of real water and simulation results for different water models.The numbers are obtain from Reference 132.while still permitting significant water transport. It is particularly true for water models suchas the TIP4P/2005 model, which have highly ordered structures at lower temperatures butless ordered at higher temperatures.We have demonstrated that the simulated water and ion transport rates through vari-ous CNTs can exhibit strong water model dependences. This observation should be kept inmind whenever MD simulations are used to investigate water and ion flow through nanoscopicchannels. A good water model should simulate real water properties as accurately as possi-ble. Recently, Vega and Abascal132 proposed a test in which many properties of water aretaken into account to evaluate the performance of a water model. All three water models(TIP4P/2005, SPC/E, and TIP3P) present in this dissertation were subjected to their test.Experimental data of real water and simulation results for different water models are tabu-lated in Table 6.1. The TIP4P/2005 model achieves the highest score in reproducing overallwater properties, especially the self-diffusion coefficient and the shear viscosity. Moreover, the766.2. Future Workexperiments29, 60 have identified quasi-one-dimensional water structures within intermediate-size nanotubes. Our simulation results (Tables 3.2, 5.5, and Figures 3.6, 3.7, 5.6) show thatthe TIP4P/2005 model tends to form stacked ring configurations, even at room temperature.Therefore, taken together, our results suggest that it would be desirable to select a watermodel such as TIP4P/2005 in future simulation studies of water and ion transport throughnanoscopic channels.6.2 Future WorkThe work in Chapter 5 generates a question: why do continuum hydrodynamic predictionsof ion entry rates agree well with the TIP3P model but fail with the TIP4P/2005 model?We proposed that the water structure near the orifice of the nanotube could be a potentialexplanation; however, this has not been justified by simulations. Studying this question mightgive further insight into entrance effect on ion conduction.As described in Chapter 1.2, the water flow rates through CNTs are observed to be fasterthan those predicted by the classical hydrodynamic equation of Poiseuille flow (Equation(1.1)). A recent paper71 summarized the results of water flow enhancement through nanotubesfrom various simulation studies. It is astonishing that these results differ by one to fiveorders of magnitude. The present dissertation at least provides some plausible explanations,including the liquid diffusivity in bulk, and the confined molecular structure of the water. Infuture work one could conduct MD simulations using appropriate water models, for example,the TIP4P/2005 model as we suggest, so as to improve the accuracy of simulation resultsthrough CNTs. Additionally, it would be interesting to investigate whether the water modeldependence exists for other types of nanoscopic channels, such as boron-nitride nanotubes(BNNTs) and aquaporins.Being effective water conduits, CNTs are regarded as a promising candidate for nano-filtration. Many studies88, 89, 107–111 have made significant efforts on modifying nanotubeproperties to make them more efficient in excluding or selecting certain ions. Because the776.2. Future Workwater structure inside nanotubes makes a substantial difference in water and ion flow rates,our results suggest that inducing the formation of polygonal ring configurations might bean important attribute to consider, which could be achieved by designing or even synthesiz-ing nanoscopic channels with desired diameters. We predict that stacked ring-like structureswithin channels should improve the separation efficiency of water from ions.Water and ion conduction through nanoscopic channels of different geometric shape is asubject of board interest. Intriguing examples are channels with the hourglass shape,156, 161, 162which have larger orifices than simple cylindrical nanotubes. Given the fact that entranceeffects are important, as demonstrated in the current dissertation, one might expect that theoverall resistance will decrease because water and ions can more easily enter the channel,however, the potential influence of structural differences between the orifices and the centralnanotube is unclear. Slablike channels are another interesting example.163 Algara-Siller etal.163 experimentally observed a “square ice” phase for water confined between two graphenesheets, and further confirmed their observations using MD simulation results. Consideringthat ordered structures within nanotubes can impede the flow of water and ions, it wouldbe interesting to learn whether this conclusion can be extrapolated to the case of slablikechannels.78Bibliography[1] P. Agre, L. S. King, M. Yasui, W. B. Guggino, O. P. Ottersen, Y. Fujiyoshi, A. Engel,and S. Nielsen, J. Physiol., 542, 3 (2002).[2] P. Agre, Angew. Chem. Int. Ed., 43, 4278 (2004).[3] T. Zeuthen, J. Membr. Biol., 234, 57 (2010).[4] A. S. Verkman, J. Cell Sci., 124, 2107 (2011).[5] T. Gonen and T. Walz, Quart. Rev. Biophys., 39, 4 (2006).[6] C. Zhao, H. Shao, and L. Chu, Colloids Surf. B, 62, 163 (2008).[7] D. F. Savage, J. D. O’Connell III, L. J. W. Miercke, J. Finer-Moore, and R. M. Stroud,Proc. Natl. Acad. Sci. U. S. A., 107, 17164 (2010).[8] Y. Wang and E. Tajkhorshid, J. Nutr., 137, 1509S (2007).[9] A. J. Yool and E. M. Campbell, Mol. Aspects Med., 33, 553 (2012).[10] B. Roux, Annu. Rev. Biophys. Biomol. Struct., 34, 153 (2005).[11] B. Corry and M. Thomas, J. Am. Chem. Soc., 134, 1840 (2012).[12] S. B. Long, E. B. Campbell, and R. MacKinnon, Science, 309, 897 (2005).[13] H. Berkefeld, B. Fakler, and U. Schulte, Physiol. Rev., 90, 1437 (2010).[14] C. Tang, Y. Zhao, R. Wang, C. He´lix-Nielsen, and A. G. Fane, Desalination, 308, 34(2013).79Bibliography[15] T. J. Jentsch, V. Stein, F. Weinreich, and A. A. Zdebik, Physiol. Rev., 82, 503 (2002).[16] S. Nielsen, J. Frøkiær, D. Marples, T. Kwon, P. Agre, and M. A. Knepper, Physiol.Rev., 82, 205 (2002).[17] A. L. George Jr., J. Clin. Invest., 115, 1990 (2005).[18] J. J. Kasianowicz, Chem. Rev., 112, 6215 (2012).[19] A. S. Verkman, Annu. Rev. Med., 63, 303 (2012).[20] A. S. Verkman, M. O. Anderson, and M. C. Papadopoulos, Nat. Rev. Drug Discov., 13,259 (2014).[21] P. R. Bandaru, J. Nanosci. Nanotechnol., 7, 1239 (2007).[22] K. Liew, Z. Lei, and L. Zhang, Compos. Struct., 120, 90 (2015).[23] Z. Han and A. Fina, Prog. Polym. Sci., 36, 914 (2011).[24] J. Prasek, J. Drbohlavova, J. Chomoucka, J. Hubalek, O. Jasek, V. Adam, and R. Kizek,J. Mater. Chem., 21, 15872 (2011).[25] A. I. Kolesnikov, J. Zanotti, C. Loong, P. Thiyagarajan, A. P. Moravsky, R. O. Loutfy,and C. J. Burnham, Phys. Rev. Lett., 93, 035503 (2004).[26] N. Naguib, H. Ye, Y. Gogotsi, A. G. Yazicioglu, C. M. Megaridis, and M. Yoshimura,Nano Lett., 4, 2237 (2004).[27] Q. Chen, J. L. Herberg, G. Mogilevsky, H. Wang, M. Stadermann, J. K. Holt, Y. Wu,Nano Lett., 8, 1902 (2008).[28] A. Das, S. Jayanthi, H. S. M. V. Deepak, K. V. Ramanathan, A. Kumar, C. Dasgupta,and A. K. Sood, ACS Nano, 4, 1687 (2010).80Bibliography[29] H. Kyakuno, K. Matsuda, H. Yahiro, Y. Inami, T. Fukuoka, Y. Miyata, K. Yanagi, Y.Maniwa, H. Kataura, T. Saito, M. Yumura, and S. Iijima, J. Chem. Phys., 134, 244501(2011).[30] E. Paineau, P. Albouy, S. Rouzie`re, A. Orecchini, S. Rols, and P. Launois, Nano Lett.,13, 1751 (2013).[31] S. D. Bernardina, E. Paineau, J.-B. Brubach, P. Judeinstein, S. Rouzie`re, P. Launois,P. Roy, J. Am. Chem. Soc., 138, 10437 (2016).[32] X. Pan, Z. Fan, W. Chen, Y. Ding, H. Luo, and X. Bao, Nat. Mater., 6, 507 (2007).[33] F. Du, L. Qu, Z. Xia, L. Feng, and L. Dai, Langmuir, 27, 8437 (2011).[34] H. Zhang, X. Pan, X. Han, X. Liu, X. Wang, W. Shen, and X. Bao, Chem. Sci., 4, 1075(2013).[35] J. Geng, K. Kim, J. Zhang, A. Escalada, R. Tunuguntla, L. R. Comolli, F. I. Allen, A.V. Shnyrova, K. Cho, D. Munoz, Y. M. Wang, C. P. Grigoropoulos, C. M. Ajo-Franklin,V. A. Frolov, and A. Noy, Nature, 514, 612 (2014).[36] W. Yu, P. Hou, F. Li, and C. Liu, J. Mater. Chem., 22, 13756 (2012).[37] J. K. Holt, H. G. Park, Y. Wang, M. Stadermann, A. B. Artyukhin, C. P. Grigoropoulos,A. Noy, and O. Bakajin, Science, 312, 1034 (2006).[38] X. Qin, Q. Yuan, Y. Zhao, S. Xie, and Z. Liu, Nano Lett., 11, 2173 (2011).[39] H. Lee, H. Kim, Y. Cho, and H. Park, Small, 10, 2653 (2014).[40] P. Pang, J. He, J. Park, P. S. Krstic´, and S. Lindsay, ACS Nano, 5, 7277 (2011).[41] W. Choi, C. L, M. Ham, S. Shimizu, and M. S. Strano, J. Am. Chem. Soc., 133, 203(2011).81Bibliography[42] W. Choi, Z. W. Ulissi, S. F. E. Shimizu, D. O. Bellisario, M. D. Ellison, and M. S.Strabo, Nat. Commun., 4, 2397 (2013).[43] L. Liu, C. Yang, K. Zhao, J. Li, and H. Wu, Nat. Commum., 4, 2989 (2013).[44] J. Wu, K. Gerstandt, H. Zhang, J. Liu, and B. J. Hinds, Nat. Nanotechnol., 7, 133(2012).[45] F. Fornasiero, H. Park, J. K. Holt, M. Stadermann, C. P. Grigoropoulos, A. Noy, andO. Bakajin, Proc. Natl. Acad. Sci. U. S. A., 105, 17250 (2008).[46] F. Taherian, V. Marcon, and N. F. A. van der Vegt, Langmuir, 29, 1457 (2013).[47] G. Hummer, J. C. Rasaiah, and J. P. Noworyta, Nature, 414, 188 (2001).[48] H. Kumar, B. Mukherjee, S. Lin, C. Dasgupta, A. K. Sood, and P. K. Maiti, J. Chem.Phys., 134, 124105 (2011).[49] H. Kumar, C. Dasgupta, and P. K. Maiti, Mol. Simul. 41, 504 (2015).[50] T. A. Pascal, W. A. Goddard, and Y. Jung, Proc. Natl. Acad. Sci. U. S. A., 108, 11794(2011).[51] S. Vaitheeswaran, J. C. Rasaiah, and G. Hummer, J. Chem. Phys., 121, 7955 (2004).[52] J. A. Garate, T. Perez-Acle, and C. Oostenbrink, Phys. Chem. Chem. Phys., 16, 5119(2014).[53] A. Waghe, J. C. Rasaiah, and G. Hummer, J. Chem. Phys., 137, 044709 (2012).[54] A. Waghe, J. C. Rasaiah, and G. Hummer, J. Chem. Phys., 117, 10789 (2002).[55] A. Striolo, Nano Lett., 6, 633 (2006).[56] H. Fang, R. Wan, X. Gong, H. Lu, and S. Li, J. Phys. D: Appl. Phys., 41, 103002(2008).82Bibliography[57] J. Ko¨finger, G. Hummer, and C. Dellago, Phys. Chem. Chem. Phys., 13, 15403 (2011).[58] K. Koga, R. D. Parra, H. Tanaka, X. C. Zeng, J. Chem. Phys., 113, 5037 (2000).[59] K. Koga, G. Gao, H. Tamaka, and X. Zeng, Nature, 412, 802 (2001).[60] Y. Maniwa, H. Kataura, M. Abe, A. Udaka, S. Suzuki, Y. Achiba, H. Kira, K. Matsuda,H. Kadowaki, and Y. Okabe, Chem. Phys. Lett., 401, 534 (2005).[61] A. Striolo, A. A. Chialvo, K. E. Gubbins, and P. T. Cummings, J. Chem. Phys., 122,234712 (2005).[62] D. Takaiwa, I. Hatano, K. Koga, and H. Tanaka, Proc. Natl. Acad. Sci. U. S. A., 105,39 (2008).[63] J. A. Thomas and A. J. H. McGaughey, Phys. Rev. Lett., 102, 184502 (2009).[64] T. Ohba, Angew. Chem. Int. Ed., 53, 8032 (2014).[65] Y. Liu, Q. Wang, T. Wu, and L. Zhang, J. Chem. Phys., 123, 234701 (2005).[66] S. Li and B. Schmidt, Phys. Chem. Chem. Phys., 17, 7303 (2015).[67] M. C. Gordillo and J. Mart´ı, Chem. Phys. Lett., 329, 341 (2000).[68] I. Hanasaki and A. Nakatani, J. Chem. Phys., 124, 144708 (2006).[69] I. Hanasaki and A. Nakatani, J. Chem. Phys., 124, 174714 (2006).[70] Q. Chen, Q. Wang, Y. Liu, and T. Wu, J. Chem. Phys., 140, 214507 (2014).[71] S. K. Kannam, B. D. Todd, J. S. Hansen, and P. J. Daivis, J. Chem. Phys., 138, 094701(2013).[72] K. Falk, F. Sedlmeler, L. Joly, R. R. Netz, and L. Bocquet, Nano Lett., 10, 4067 (2010).[73] S. Joseph and N. R. Aluru, Nano Lett., 8, 452 (2008).83Bibliography[74] C. Won, S. Joseph, and N. R. Aluru, J. Chem. Phys., 125, 114701 (2006).[75] P. Sahu and S. M. Ali, J. Chem. Phys., 143, 184503 (2015).[76] D. Lu, Phys. Chem. Chem. Phys., 15, 14447 (2013).[77] J. Li, X. Gong, H. Lu, D. Li, H. Fang, and R. Zhou, Proc. Natl. Acad. Sci. U. S. A.,104, 3687 (2007).[78] X. Gong, J. Li, H. Lu, R. Wan, J. Li, J. Hu, and H. Fang, Nat. Nanotechnol., 2, 709(2007).[79] L. Liu, Y. Qiao, and X. Chen, Appl. Phys. Lett., 92, 101927 (2008).[80] X. Li, Y. Shi, Y. Yang, H. Du, R. Zhou, and Y. Zhao, J. Chem. Phys., 136, 175101(2012).[81] S. Hou, Z. Shen, X. Zhao, and Z. Xue, Chem. Phys. Lett., 373, 308 (2003).[82] D. Lu, Y. Li, U. Ravaioli, and K. Schulten, J. Phys. Chem. B, 109, 11461 (2005).[83] R. Wan, J. Li, H. Lu, and H. Fang, J. Am. Chem. Soc., 127, 7166 (2005).[84] L. Wu, F. Wu, J. Kou, and H. Lu, Phys. Rev. E, 83, 061913 (2011).[85] M. Melillo, F. Zhu, M. A. Snyder, and J. Mittal, J. Phys. Chem. Lett., 2, 2978 (2011).[86] I. Moskowitz, M. A. Snyder, and J. Mittal, J. Chem. Phys., 141, 18C532 (2014).[87] F. Ranazani and F. Ebrahimi, J. Phys. Chem. C, 120, 12871 (2016).[88] B. Corry, Energy Environ. Sci., 4, 751 (2011).[89] Z. E. Hughes, C. J. Shearer, J. Shapter, and J. D. Gale, J. Phys. Chem. C, 116, 24943(2012).[90] C. Peter and G. Hummer, Biophys. J., 89, 2222 (2005).84Bibliography[91] L. A. Richards, A. I. Scha¨fer, B. S. Richards, and B. Corry, Small, 8, 1701 (2012).[92] H. Liu, S. Murad, and C. J. Jameson, J. Chem. Phys., 125, 084713 (2006).[93] Q. Shao, L. Huang, J. Zhou, L. Lu, L. Zhang, X. Lu, S. Jiang, K. E. Gubbins, and W.Shen, Phys. Chem. Chem. Phys., 10, 1896 (2008).[94] E. I. Calixte, O. N. Samoylova, and K. L. Shuford, Phys. Chem. Chem. Phys., 18, 12204(2016).[95] B. Corry, J. Phys. Chem. B, 112, 1427 (2008).[96] O. N. Samoylova, E. I. Calixte, and K. L. Shuford, J. Phys. Chem. C, 119, 1659 (2015).[97] T. A. Beu, J. Chem. Phys., 132, 164513 (2010).[98] Z. He, B. Corry, X. Lu, and J. Zhou, Nanoscale, 6, 3686 (2014).[99] J. Liu, G. Shi, P. Guo, J. Yang, and H. Fang, Phys. Rev. Lett., 115, 164502 (2015).[100] M. C. F. Wander and K. L. Shuford, J. Phys. Chem. C, 114, 20539 (2010).[101] R. Qiao and N. R. Aluru, Nano Lett., 3, 1013 (2003).[102] T. Sumikama, S. Saito, and I. Ohmine, J. Phys. Chem. B, 110, 20671 (2006).[103] T. Sumikama, S. Saito, and I. Ohmine, J. Chem. Phys., 139, 165106 (2013).[104] J. Su and D. Huang, J. Phys. Chem. C, 120, 11245 (2016).[105] T. A. Beu, J. Chem. Phys., 135, 044516 (2011).[106] J. Azamat, A. Khataee, and S. Joo, RSC Adv. 5, 25097 (2015).[107] C. Song and B. Corry, J. Phys. Chem. B, 113, 7642 (2009).[108] H. Liu, C. J. Jameson, and S. Murad, Mol. Simul., 34, 169 (2008).[109] Q. Shao, J. Zhou, L. Lu, X. Lu, Y. Zhu, and S. Jiang, Nano Lett., 9, 989 (2009).85Bibliography[110] D. Tang and D. Kim, Chem. Phys., 428, 14 (2014).[111] X. Gong, J. Li, K. Xu, J. Wang, and H. Yang, J. Am. Chem. Soc., 132, 1873 (2010).[112] J. D. Bernal and R. H. Fowler, J. Phys. Chem., 1, 515 (1933).[113] W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey, and M. L. Klein, J.Chem. Phys. 79, 926 (1983).[114] H. J. C. Berendsen, J. R. Grigera, and T. P. Straatsma, J. Phys. Chem. 91, 6269 (1987).[115] P. Ren and J. W. Ponder, J. Phys. Chem. B, 107, 5933 (2003).[116] J. L. F. Abascal and C. Vega, J. Chem. Phys. 123, 234505 (2005).[117] W. Yu, P. E. M. Lopes, B. Roux, and A. D. MacKerell, Jr., J. Chem. Phys., 138, 034508(2013).[118] B. Guillot, J. Mol. Liq., 101, 219 (2002).[119] S. Tanizaki, J. Mavri, H. Partridge, and P. C. Jordan, Chem. Phys., 246, 37 (1999).[120] S. Zhu and C. F. Wong, J. Chem. Phys., 98, 8892 (1993).[121] D. M. Ferguson, J. Comput. Chem., 16, 501 (1995).[122] Y. Wu, H. L. Tepper, and G. A. Voth, J. Chem. Phys., 124, 024503 (2006).[123] M. A. Gonza´lez and J. L. F. Abascal, J. Chem. Phys., 135, 224516 (2011).[124] J. Wang, P. Cieplak, Q. Cai, M. Hsieh, J. Wang, Y. Duan, and R. Luo, J. Phys. Chem.B, 116, 7999 (2012).[125] A. Alexiadis and S. Kassinos, Chem. Rev., 108, 5014 (2008).[126] A. Alexiadis and S. Kassinos, Chem. Eng. Sci., 63, 2793 (2008).[127] A. Alexiadis and S. Kassinos, Mol. Simul., 34, 671 (2008).86Bibliography[128] H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, and J. Hermans, Inter-molecular Forces (Reidel, Dordrecht, 1981), p. 331.[129] Y. Nakamura and T. Ohno, Mater. Chem. Phys., 132, 682 (2012).[130] S. W. Rick, J. Chem. Phys., 120, 6085 (2004).[131] H. Kumar, C. Dasgupta, and P. K. Maiti, RSC Adv., 5, 1893 (2015).[132] C. Vega and J. L. F. Abascal, Phys. Chem. Chem. Phys., 13, 19663 (2011).[133] S. Andreev, D. R. Reichman, and G. Hummer, J. Chem. Phys., 123, 194502 (2005).[134] T. Darden, D. York, and L. Pedersen, J. Chem. Phys., 98, 10089 (1993).[135] U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee, and L. G. Pedersen, J.Chem. Phys., 103, 8577 (1995).[136] Y. Duan, C. Wu, S. Chowdhury, M. C. Lee, G. Xiong, W. Zhang, R. Yang, P. Cieplak,R. Luo, T. Lee, J. Caldwell, J. Wang, and P. Kollman, J. Comput. Chem. 24, 1999(2003).[137] I. Joung and T. E. Cheatham, III, J. Phys. Chem. B, 112, 9020 (2008).[138] D. van der Spoel, E. Lindahl, B. Hess, G. Groenhof, A. E. Mark, and H. J. C. Berendsen,J. Comput. Chem., 26, 1701 (2005).[139] G. Bussi, D. Donadio, and M. Parrinello, J. Chem. Phys., 126, 014101 (2007).[140] F. Zhu, E. Tajkhorshid, and K. Schulten, Biophys. J., 86, 50 (2004).[141] F. Zhu, E. Tajkhorshid, and K. Schulten, Biophys. J., 83, 154 (2002).[142] M. E. Suk and N. R. Aluru, J. Phys. Chem. Lett., 1, 1590 (2010).[143] L. Wang, R. S. Dumont, and J. M. Dickson, J. Chem. Phys., 137, 044102 (2012).[144] A. Luzar and D. Chandler, J. Chem. Phys., 98, 8160 (1993).87Bibliography[145] G. M. Torrie and J. P. Valleau, J. Comput. Phys., 23, 187 (1977).[146] J. Ka¨stner, WIREs Comput. Mol. Sci., 1, 932 (2011).[147] B. Roux, Comput. Phys. Comm., 91, 275 (1995).[148] L. Liu and G. N. Patey, J. Chem. Phys., 141, 18C518 (2014).[149] L. Liu and G. N. Patey, J. Chem. Phys., 144, 184502 (2016).[150] T. A. Hilder, D. Gordon, and S. Chung, Small, 5, 2183 (2009).[151] J. Goldsmith and C. C. Martens, J. Phys. Chem. Lett., 1, 528 (2010).[152] R. M. Lynden-Bell and J. C. Rasaiah, J. Chem. Phys., 105, 9266 (1996).[153] G. Lakatos and G. N. Patey, J. Chem. Phys., 126, 024703 (2007).[154] T. B. Sisan and S. Lichter, Microfluid. Nanofluid., 11, 787 (2011).[155] J. H. Walther, K. Ritos, E. R. Cruz-Chu, C. M. Megaridis, and P. Koumoutsakos, NanoLett., 13, 1910 (2013).[156] S. Gravelle, L. Joly, C. Ybert, and L. Bocquet, J. Chem. Phys., 141, 18C526 (2014).[157] R. A. Sampson, Phil. Trans. R. Soc. Lond. A, 182, 449 (1891).[158] L. Liu and G. N. Patey, J. Chem. Phys., 146, 074502 (2017).[159] M. Thomas and B. Corry, Phil. Trans. R. Soc. A, 374, 20150220 (2016).[160] R. Mancinelli, A. Botti, F. Bruni, M. A. Ricci, and A. K. Soper, J. Phys. Chem. B,111, 13570 (2007).[161] D. Tang, L. Li, M. Shahbabaei, Y. Yoo, and D. Kim, Materials, 8, 7257 (2015).[162] M. Shahbabaei and D. Kim, Chem. Phys., 477, 24 (2016).88[163] G. Algara-Siller, O. Lehtinen, F. Wang, R. R. Nair, U. Kaiser, H. Wu, A. K. Geim, andI. V. Grigorieva, Nature, 519, 443 (2015).[164] P. Ewald, Ann. Phys., 369, 253 (1921).[165] A. Y. Toukmaji and J. A. Board Jr., Comput. Phys. Comm., 95, 73 (1996).[166] R. W. Hockney, S. P. Goel, J. W. Eastwood, J. Comput. Phys., 14, 148 (1974).[167] B. Hess, H. Bekker, H. J. C. Berendsen, J. G. E. M. Fraaije, J. Comput. Chem., 18,1463 (1997).[168] J. G. Kirkwood, J. Chem. Phys., 3, 300 (1935).[169] A. Kalra, S. Garde, G. Hummer, Proc. Natl. Acad. Sci. U. S. A., 100, 10175 (2003).[170] W. D. Nicholls, M. K. Borg, D. A. Lockerby, and J. M. Reese, Microfluid Nanofluid,12, 257 (2012).[171] M. Thomas and B. Corry, Microfluid. Nanofluid., 18, 41 (2015).[172] H. J. C. Berendsen, J. P. M. Postma, W. F. Vangunsteren, A. Dinola, and J. R. Haak,J. Chem. Phys., 81, 3684 (1984).89Appendix AMolecular DynamicsA.1 Ewald SummationAs described in Chapter 2.3, the pair electrostatic potential, uC,ij, falls off slowly, such that thecontributions of long-range interactions to the total potential energy are not negligible. Withthe application of the periodic boundary conditions, the total Coulombic potential energy ofpoint charges can be expressed asUC,total =18πǫ0∑i∑j∑n′qiqj|rij + nL| , (A.1)where rij is the vector from site i to j, n is a vector (nx, ny, nz) specifying the translationalindex of the image cell, L is a set of coefficients defining the dimensions of the simulation cell(x, y, z), and the prime indicates that the i = j terms for n = (0, 0, 0) are omitted to avoidcalculating self interactions. However, this direct sum is not only slowly but also conditionallyconvergent.Ewald summation164, 165 is introduced to calculate electrostatic interactions. Unlike thedirect summation method, which uses discrete Dirac delta functions to describe the chargedistributions, the Ewald summation method adds a Gaussian charge distribution of differentsign to the original charge distributions, and then subtracts the same Gaussian distribution,with a form ofρi(r) =qiα3π32e−α2r2 , (A.2)where r is the position relative to the central point charge, and α, which is called the Ewald90A.1. Ewald Summationconvergence parameter, is a positive parameter that determines the width of the distribution.Thereafter, the conditionally convergent sum in Equation (A.1), is reconstructed into threecomponents, two quickly convergent series, and a constant termUC,total = Ureal + Urec + Uconst . (A.3)The two series, Ureal and Urec, are calculated in real spaceUreal =18πǫ0∑i∑j∑n′qiqjerfc(α|rij + nL|)|rij + nL| , (A.4)and reciprocal spaceUrec =12πV∑i∑jqiqj∑n′exp(− (πn/α)2 + 2πinrij)n2, (A.5)respectively, where “erfc” indicates the complementary error function. Note that in Equation(A.4), a large value of αmakes the real space series converge more rapidly, whereas in Equation(A.5), a small value of α makes the reciprocal space series converge more rapidly. It isimportant that the simulation cell should be exactly neutral. Otherwise, the series will notbe convergent. The constant term Uconst is a correction termUconst = − α4πǫ0√π∑i∑αq2iα , (A.6)which cancels out the interactions of each of introduced artificial charge distributions withitself. The computer time required for the sum of the reciprocal component in the traditionalEwald method increases as N2, where N is the number of point charges in a simulation cell,therefore, it is not realistic for large systems.In GROMACS, the particle mesh Ewald (PME) method,134, 135 which is based on fastFourier transform (FFT), is used to accelerate the summation of the reciprocal part. Inthis method, Ureal in Equation (A.4) is treated as in the traditional Ewald method. Urecin Equation (A.5) is handled by a particle-mesh technique, which converts point charges91A.2. System Evolutioninto a mesh of density values. The electrostatic potentials are calculated by solving Poissonequations using FFT techniques from the mesh field, and so are the forces on point charges.The computer time required for the PME method increases as Nlog(N), which is much fasterthan the traditional Ewald method when N is large.A.2 System EvolutionIn our simulations, we employ the leap-frog algorithm166 as the MD integrator. This algorithmuses atom coordinates r at time t and atom velocities v at time t−δt/2. The method updatesvelocities and coordinates using the equations of motionv(t+δt2) = v(t− δt2) + δtM−1f(t) , (A.7)r(t+ δt) = r(t) + δtv(t+δt2) , (A.8)where M is a diagonal matrix containing the masses of the particles, and f(t) is the forcevector, which is given byf(t) = −∇(uLJ(t) + uC(t)) . (A.9)With rigid models, the distances between pairs of atoms which form chemical bonds areconstrained to definite bond lengths d1, d2, . . . di. Therefore, a linear constraint solver(LINCS) algorithm167 is implemented in GROMACS. This algorithm corrects unconstrainedvelocities and coordinates using a two-step method. First, taking constraints into considera-tion, the constrained equations of motion arev(t+δt2) =(I − T (t)B(t))(v(t− δt2) + δtM−1f(t))− T (t)(B(t)r(t)− d)δt, (A.10)r′(t+ δt) =(I − T (t)B(t))(r(t) + δtv(t− δt2) + δt2M−1f(t))+ T (t)d , (A.11)92A.3. Thermostatwhere I is the identity matrix, B(t) is a matrix containing the directions of the constraints,d is a vector containing the prescribed bond lengths, and T (t) is a matrix defined by theexpression M−1B(t)T(B(t)M−1B(t)T)−1. However, Equation (A.11) sets the projectionsof new bonds onto the directions of old bonds to the prescribed lengths, instead of setting thelengths of new bonds l1, l2, . . . li themselves. Therefore, a second stepr(t+ δt) =(I − T (t)B(t))r′(t+ δt) + T (t)p(t) , (A.12)is required to update coordinates, where p(t) is a projection vector for the correction. Theelements of p(t) are given bypi =√2d2i − l2i . (A.13)A.3 ThermostatIn this dissertation, we used a refined velocity rescaling algorithm proposed by Bussi et al.139for sampling the canonical distribution. It has a general procedure as follows:(1) Evolve the system for a single time step, as described in A.2.(2) Calculate the system kinetic energy.(3) Evolve the kinetic energy for a single time step using an auxiliary stochastic dynamics.The auxiliary dynamics obeys the expressiondK = (K0 −K) dtτ+ 2√KK0NfdW√τ, (A.14)where K is the system kinetic energy, K0 is the desired kinetic energy, dt is the time step, τ isan arbitrary parameter with the dimension of time, Nf is the number of degrees of freedom,and dW is a stochastic Wiener process (a random number). The parameter τ , sometimescalled the relaxation time, defines the intensity of temperature coupling. A larger τ means ittakes longer to achieve the given kinetic energy. In our simulations, a τ of 0.1 ps was used.93A.4. Potential of Mean Force(4) Rescale velocities in the manner of(v(t+δt2))2=(v′(t+δt2))2+ 2M−1dK , (A.15)so as to update the kinetic energy, where v′(t + δt2 ) and v(t +δt2 ) are velocities before andafter rescaling.A.4 Potential of Mean ForceThe potential of mean force (PMF) of a system with N particles is the potential, if a set ofparticles 1 . . .n keeps a fixed configuration, that gives the average force over all configurationsof all remaining n+ 1 . . .N particles acting on a particle i,168−∇iw =∫(−∇jU)e−U(r)kBT drn+1 . . . drN∫e−U(r)kBT drn+1 . . . drN, (A.16)where w is the PMF (−∇iw is the average force), U is the system potential energy, kB is theBoltzmann constant, T is the absolute temperature, and r is the particle coordinates in thephase space.The PMF w(η) as a function of the reaction coordinate η is defined from the averagedistribution functionw(η) = w(η∗)− kBT ln[ 〈ρ(η)〉〈ρ(η∗)〉], (A.17)where η∗ and w(η∗) are arbitrary constants. The average distribution function along thecoordinate η, 〈ρ(η)〉, is obtained from a Boltzmann weighted average〈ρ(η)〉 =∫δ(η′[r]− η)e−U(r)kBT dr∫e−U(r)kBT dr, (A.18)where η′[r] is a function which allows the disturbance of a few degrees of freedom, and94A.4. Potential of Mean Forceδ(η′[r]− η) is the Dirac delta function for the coordinate η.The measure of a thermodynamic property obtained from average distribution functions isaccurate provided that the system is ergodic. However, a high energy barrier along the reactioncoordinate may cause poor sampling in the high energy parts of phase space, therefore, itis unrealistic to compute w(η) or 〈ρ(η)〉 directly from straight MD simulations. Umbrellasampling145, 146 is a technique to overcome the sampling problem by introducing artificial biaspotentials, Ui(η). In practice multiple simulations (windows) are conducted with differentbiased potentials Ubias(r), denoted by the subscript i. The unbiased PMF can be readilyevaluated by biased average distribution functionsw(η)i = w(η∗)− kBT ln[〈ρ(η)〉biasi〈ρ(η∗)〉]− Ubias(r) + Fi , (A.19)where Ubias(r) is the biased potentialUbias(r) = U(r) + Ui(η) , (A.20)and Fi is a constant associated with Ui(η) bye−FikBT =〈e−Ui(η)kBT〉. (A.21)After obtaining individual w(η)i the free energy curves are recovered with the weightedhistogram analysis method (WHAM).147 The global 〈ρ(η)〉 can be written as〈ρ(η)〉 =∑ini 〈ρ(η)〉biasi∑inie−wi(η)−FikBT, (A.22)where n is the number of independent data points of individual simulations. The constantsFi required in the previous equation are determined frome−FikBT =∫〈ρ(η)〉 e−wi(η)kBT dη . (A.23)95A.4. Potential of Mean Forcewhich again needs the input of 〈ρ(η)〉. Practically, Equations (A.22) and (A.23) are solvedthrough an iteration procedure so that 〈ρ(η)〉 and the set of Fi are self-consistent.96Appendix BDiscussion of Some SimulationDetailsB.1 Cell DimensionTo verify that the dimensions of our simulation cell are not seriously influencing the results,additional simulations were performed with higher cell reservoirs (4.004 versus 2.004 nm), andthe results obtained are given in Table B.1. We note that the results agree within standarddeviations. It also has been reported95, 142, 169, 170 that the length of CNT has a negligibleeffect on the water conduction rate.B.2 ThermostatAs demonstrated by Thomas and Corry,171 the thermostat employed can influence waterflow dynamics through nanoscopic channels. They found that temperature control algorithmsthat introduce random forces, for example, the Langevin method, can lead to flow rates whichdepend on the length of the CNT, whereas velocity rescaling algorithms do not. Although itis not clear which type of thermostat more closely resembles a real experimental situation,we used the velocity rescaling algorithm presented by Bussi et al.139 We also examined theBerendsen thermostat,172 which is one of the velocity rescaling thermostats investigated byThomas and Corry,171 to test for any influence on the flow rate. The results are tabulatedin Table B.2 and suggest that the flow rate is not unduly sensitive to details of the velocityrescaling thermostat employed.97B.3. Implementation of Pressure DifferenceB.3 Implementation of Pressure DifferenceIn Chapter 2.5 we discussed that an improper selection of the region where the external forcefex is applied could adversely influence the fluid dynamics. To validate our choice of regionthickness, which is 0.2 nm, we further tested two other thickness values: 0.4 nm and 2.0 nm(full reservoir). The results are summarized in Table B.3, and confirm that our applicationis reasonable, while the full reservoir method of Reference 141 obviously alters the fluiddynamics.To work more efficiently with the GROMACS package, we arbitrarily defined the timeinterval of subset update as 10 ps. Another two periods, 5 ps and 25 ps, were investigatedunder the same conditions. The results are displayed in Table B.4, demonstrating that theselection of the time interval does not have a discernible effect on water flow rates.We met an unexpected issue for simulations involving ions. If we applied the externalforce to both water molecule and ions within the selected region, it resulted in unequal NaClconcentrations at steady state, with the feed reservoir having a higher concentration thanthe downstream reservoir. For example, a feed solution initially at 1 mol L−1 ends up atapproximately 1.6 mol L−1 at steady state. Also, the observed ion flow rates increase, as shownin Table B.5. We believe that this is a systematic artifact of periodic boundary conditionsand finite size. Concentration separation did not occur if the force was applied only to watermolecules. Interestingly, when we carried out simulations for a 1.6 mol L−1 solution withoutthe exerting force on the ions, we obtained ion flow rates similar to those for a 1 mol L−1solution applying force on the ions (Table B.5). It suggests that the flow rate increase is mainlydue to the ion concentration increase in the feed reservoir, rather than to some fundamentalchanges in flow dynamics.98B.3. Implementation of Pressure DifferenceWater model Rflow,2.0nm (ns−1) Rflow,4.0nm (ns−1)TIP4P/2005 22 (3) 23 (2)SPC/E 34 (2) 30 (3)TIP3P 59 (2) 61 (4)Table B.1: Effect of the cell dimension on water flow rates. The simulations were conductedat 300 K for a (6,6) CNT. Rflow is the average water flow rate. The numbers in parenthesisare estimated standard deviations.Water model Rflow,Bussi (ns−1) Rflow,Berendsen (ns−1)TIP4P/2005 22 (3) 20 (3)SPC/E 34 (2) 31 (4)TIP3P 59 (2) 61 (2)Table B.2: Effect of the thermostat employed on water flow rates. The simulations wereconducted at 300 K for a (6,6) CNT. Rflow is the average water flow rate. The numbers inparenthesis are estimated standard deviations.Water model Rflow,0.2nm (ns−1) Rflow,0.4nm (ns−1) Rflow,2.0nm (ns−1)TIP4P/2005 22 (3) 22 (2) 48 (2)SPC/E 34 (2) 32 (2) 68 (1)TIP3P 59 (2) 63 (1) 131 (5)Table B.3: Effect of the thickness of the force-exerted region on water flow rates. The simu-lations were conducted at 300 K for a (6,6) CNT. Rflow is the average water flow rate. Thenumbers in parenthesis are estimated standard deviations.Water model Rflow,10ps (ns−1) Rflow,5ps (ns−1) Rflow,25ps (ns−1)TIP4P/2005 22 (3) 21 (1) 22 (2)SPC/E 34 (2) 35 (1) 36 (4)TIP3P 59 (2) 64 (3) 59 (1)Table B.4: Effect of the update frequency of the force exerted region on water flow rates. Thesimulations were conducted at 300 K for a (6,6) CNT. Rflow is the average water flow rate.The numbers in parenthesis are estimated standard deviations.99B.3. Implementation of Pressure DifferenceForce onionscinitial(mol L−1)csteady(mol L−1)Water modelRflow (ns−1)Na+ Cl−No 1 1TIP4P/2005 0.4 (0.2) 0.2 (0.2)TIP3P 3.1 (0.4) 2.2 (0.5)Yes 1 ∼1.6 TIP4P/2005 0.8 (0.4) 0.6 (0.3)TIP3P 5.2 (0.2) 4.1 (0.4)No 1.6 1.6TIP4P/2005 0.6 (0.4) 0.4 (0.2)TIP3P 4.8 (0.3) 3.7 (0.3)Table B.5: Effect of the external force on ion flow rates. cinitial and csteady are concentrationsof the solution in the feed reservoir before and after steady states were reached, respectively.Rflow is the average ion flow rate. The numbers in parenthesis are estimated standard devia-tions.100
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- A molecular dynamics investigation of water and ion...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
A molecular dynamics investigation of water and ion transport through model carbon nanotubes Liu, Lin 2017
pdf
Page Metadata
Item Metadata
Title | A molecular dynamics investigation of water and ion transport through model carbon nanotubes |
Creator |
Liu, Lin |
Publisher | University of British Columbia |
Date Issued | 2017 |
Description | In this dissertation, we investigate water and ion transport through carbon nanotubes using molecular dynamics simulations. Specifically, we examine how different water models influence the simulated conduction rates. We consider three common water models, which are TIP4P/2005, SPC/E, and TIP3P, and observe that water flow rates through the same nanotube are strikingly different amongst the different water models. Also, the water flow rate dependence on temperature fits an Arrhenius-type equation over a temperature range from 260 to 320 K. We provide evidence that there are two factors which determine the conduction rate: the bulk fluid mobility, and the molecular structure of confined water. For narrow nanotubes, for example, a (6,6) nanotube, where water only forms a single-file configuration, the first factor can largely account for the flow rate differences. In this case, we show that the conduction rate correlates with the diffusion coefficient of bulk water. Our simulation results are well described by continuum hydrodynamics as well. The factor of bulk fluid mobility is still important in the water conduction through intermediate-size nanotubes, such as a (9,9) nan- otube. Also, the formation of complex configurations within such nanotubes can impede the transport rate by influencing the mode of water conduction. The ordered structure occurring within nanotubes can also explain the differences between simulation results and continuum hydrodynamics predictions. Hence, both factors decide the water conduction rates through intermediate-size nanotubes. Moreover, we demonstrate that the ion flow rate depends on the viscosity of the bulk solution, as well as the water structure within the nanotubes, together with the ion size. In particular, at lower temperatures complex water configurations act to impede ion transport while still allowing water to flow at a significant rate. In general, our efforts on this issue are of importance for future simulation studies investigating water and ion conduction through nanoscopic channels. This dissertation might also prove useful in designing more efficient nanoscopic conduits for future experimental studies. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2017-03-15 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0343236 |
URI | http://hdl.handle.net/2429/60917 |
Degree |
Doctor of Philosophy - PhD |
Program |
Chemistry |
Affiliation |
Science, Faculty of Chemistry, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2017-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2017_may_liu_lin.pdf [ 2.27MB ]
- Metadata
- JSON: 24-1.0343236.json
- JSON-LD: 24-1.0343236-ld.json
- RDF/XML (Pretty): 24-1.0343236-rdf.xml
- RDF/JSON: 24-1.0343236-rdf.json
- Turtle: 24-1.0343236-turtle.txt
- N-Triples: 24-1.0343236-rdf-ntriples.txt
- Original Record: 24-1.0343236-source.json
- Full Text
- 24-1.0343236-fulltext.txt
- Citation
- 24-1.0343236.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0343236/manifest