UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

A molecular dynamics investigation of ice nucleation induced by electric fields Yan, Jingyi 2014

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata


24-ubc_2014_september_yan_jingyi.pdf [ 5.85MB ]
JSON: 24-1.0167534.json
JSON-LD: 24-1.0167534-ld.json
RDF/XML (Pretty): 24-1.0167534-rdf.xml
RDF/JSON: 24-1.0167534-rdf.json
Turtle: 24-1.0167534-turtle.txt
N-Triples: 24-1.0167534-rdf-ntriples.txt
Original Record: 24-1.0167534-source.json
Full Text

Full Text

A molecular dynamics investigation ofice nucleation induced by electric fieldsbyJingyi YanB.Sc., Shanghai Jiao Tong University, 2008Ph.D., The University of British Columbia, 2014A 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)July 2014c© Jingyi Yan 2014AbstractThis thesis aims to understand the influence of electric fields on ice nucle-ation. Molecular dynamics simulations are employed to investigate hetero-geneous ice nucleation induced by electric fields, and why external electricfields promote freezing in liquid water models.The first project considers heterogeneous ice nucleation in systems, wherewater molecules experience an electric field in a narrow region over an entiresurface. The specific focus is ice nucleation and growth processes. Differentwater models are considered, and the influences of temperature and fieldparameters are examined. We find no qualitative difference between the twowater models. By analyzing structure, we show that a ferroelectric cubic icelayer freezes inside the field region, and unpolarized ice grows beyond thefield region, at temperatures not far below the melting point.We explore ice nucleation by electric field bands, which act only overa portion of a surface. Field bands of different geometry nucleate ice,provided that the band is sufficiently large. Analysis of different systemsreveals that ice strongly prefers to grow at the (111) crystal plane of cubicice, and that ice nucleated by field bands usually grows as a mixture ofcubic and hexagonal ice. Our results suggest that local electric fields couldplay a major role in heterogeneous ice nucleation, particularly for roughparticles with many surface structural variations, that serve as ice nuclei inthe environment.We also investigate the electrofreezing of water subject to a uniform field.The aim is to obtain an understanding of why electric fields facilitate icenucleation. It is shown that the melting point of water increases significantlywhen water is polarized by a field. The increased melting point is mainlydue to the favourable interaction of near perfectly polarized cubic ice withiiAbstractthe applied field. Relevant to the mechanism of heterogeneous ice nucleationby local surface fields, our results suggest that local fields effectively increasethe degree of supercooling of locally polarized liquid. This decreases the sizeof the critical nucleus in the region influenced by the field, facilitating icenucleation.iiiPrefaceChapter 3 and Chapter 4 are co-authored, peer-reviewed journal articles.Chapter 5 has been submitted to a peer-reviewed journal as a co-authoredjournal article. The details of my contributions to each chapter mentionedabove are listed below.A version of Chapter 3 has been published by J. Y. Yan and G. N.Patey, Heterogeneous Ice Nucleation Induced by Electric Fields, J. Phys.Chem. Lett. 2(20), 2555, (2011); Molecular Dynamics Simulations of IceNucleation by Electric Fields, J. Phys. Chem. A. 116(26), 7057, (2012).My contributions are as follows:• Designed the research project with my supervisor.• Formulated the research questions and wrote the computation pro-grams.• Performed all of the simulations.• Performed all of the data analysis.• Prepared all of the figures in the publications.• Wrote the first draft of the publications.• Writing of the text for the publications was shared with my supervisor.A version of Chapter 4 has been published by J. Y. Yan and G. N. Patey,Ice nucleation by electric surface fields of varying range and geometry, J.Chem. Phys. 139, 144501, (2013). My contributions are as follows:• Designed the research project with my supervisor.ivPreface• Performed all of the simulations.• Performed all of the data analysis.• Prepared all of the figures in the publication.• Wrote the first draft of the publication.• Writing of the text for the publication was shared with my supervisor.A version of Chapter 5 has been submitted by J. Y. Yan, S. D. Overduin,and G. N. Patey, Understanding electrofreezing in water simulations, J.Chem. Phys. (submitted, May, 2014). My contributions are as follows:• Designed the research project with my supervisor.• Performed all the simulations with Dr. S. D. Overduin.• Performed all of the data analysis with Dr. S. D. Overduin.• Prepared the tables and figures in the publication.• Wrote the first draft of the publication.• Shared writing of the text for the publication with Dr. S. D. Overduinand my supervisor.• Additional contributions from co-authors:– Dr. S. D. Overduin did the data analysis for the tetrahedral or-der parameter distributions, structure factors, and partial radialdistribution functions.– Dr. S. D. Overduin performed the simulations in Table 5.3, andcontributed to the discussion on the influence of electric field.vTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . xvii1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Nucleation in the Atmosphere and in Biological Systems . . 11.2 Nucleation Theory . . . . . . . . . . . . . . . . . . . . . . . . 31.3 Field Induced Ice Nucleation . . . . . . . . . . . . . . . . . . 61.4 Thesis Motivation and Objectives . . . . . . . . . . . . . . . 101.5 Thesis Overview . . . . . . . . . . . . . . . . . . . . . . . . . 152 Models and Simulation Methods . . . . . . . . . . . . . . . . 162.1 Molecular Models . . . . . . . . . . . . . . . . . . . . . . . . 162.2 Interaction Potentials . . . . . . . . . . . . . . . . . . . . . . 202.3 Models of Electric Fields . . . . . . . . . . . . . . . . . . . . 242.4 System Evolution . . . . . . . . . . . . . . . . . . . . . . . . 292.4.1 The Gear Predictor-Corrector Method for TranslationalMotion . . . . . . . . . . . . . . . . . . . . . . . . . . 292.4.2 Quaternions and Rotational Motion . . . . . . . . . . 32viTable of Contents2.4.3 Simulation Constraints . . . . . . . . . . . . . . . . . 352.5 Comparison with Open Source Packages . . . . . . . . . . . . 372.6 Data Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . 392.6.1 Radial Distribution Function . . . . . . . . . . . . . . 392.6.2 The CHILL Ice Detector . . . . . . . . . . . . . . . . 402.6.3 Orientational Order Parameter . . . . . . . . . . . . . 413 Heterogeneous Ice Nucleation Induced by Surface ElectricFields . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 443.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 443.2 The Model and Simulation Method . . . . . . . . . . . . . . 463.3 Results and Discussion . . . . . . . . . . . . . . . . . . . . . 493.3.1 Ice Nucleation Induced by a Surface field . . . . . . . 493.3.2 Ice Nucleation in Systems of Varying Size and Dimen-sion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 543.3.3 Influence of Temperature . . . . . . . . . . . . . . . . 613.3.4 Influence of Field Strength and Extent . . . . . . . . 663.4 Summary and Conclusions . . . . . . . . . . . . . . . . . . . 704 Ice Nucleation by Electric Surface Fields of Varying Rangeand Geometry . . . . . . . . . . . . . . . . . . . . . . . . . . . . 744.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 744.2 Models and Methods . . . . . . . . . . . . . . . . . . . . . . 784.3 Results and Discussion . . . . . . . . . . . . . . . . . . . . . 814.3.1 Full Surface Field . . . . . . . . . . . . . . . . . . . . 814.3.2 Partial Surface Fields . . . . . . . . . . . . . . . . . . 854.4 Summary and Conclusions . . . . . . . . . . . . . . . . . . . 925 Understanding Electrofreezing in Water Simulations . . . 945.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 945.2 Simulation and Analysis Method . . . . . . . . . . . . . . . . 965.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1005.3.1 Melting Temperature . . . . . . . . . . . . . . . . . . 1005.3.2 Ice Clusters and Nucleation . . . . . . . . . . . . . . 104viiTable of Contents5.3.3 Thermodynamic Functions . . . . . . . . . . . . . . . 1115.3.4 Structural Changes with Field . . . . . . . . . . . . . 1165.4 Summary and Conclusions . . . . . . . . . . . . . . . . . . . 1196 Conclusions and Perspective . . . . . . . . . . . . . . . . . . . 1236.1 Simulation Results on Ice Nucleation and Electrofreezing . . 1236.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . 126Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128viiiList of Tables2.1 Water model and surface parameters used in the simulations.Note that ε and σ are the energy and distance parameters forthe interaction calculation, r1 is the oxygen-hydrogen bonddistance, r2 is the distance between the oxygen and the neg-ative charges (OM or OL), q is the charge on each site, andθ and φ are angles associated with different sites in watermolecules. The subscripts HOH, MOH and LOL indicate theangles between two hydrogen atoms, between the negativecharge M and the hydrogen atom, and between the two lonepair charge sites. Note that L is short for lone pair (LP). . . 183.1 The systems considered in tests of possible finite size effects.In the third column, the first and second pairs of integersgive the number of trials and the run times for the six-siteand TIP4P/Ice models, respectively. The number of trialsfor which ice nucleation was observed in the given simulationtime is given in brackets. . . . . . . . . . . . . . . . . . . . . . 553.2 Temperatures considered for the six-site and TIP4P/Ice mod-els. Yes and No indicate whether ice nucleation and growthwas observed or not at the particular temperature. . . . . . . 633.3 The dependence of ice nucleation on the parameters Emaxand c in Equation (3.1). The simulations were at 270 K forthe six-site model, and at 260 K for TIP4P/Ice. Yes and Noindicate whether ice nucleation and growth was observed ornot for the given field parameters. . . . . . . . . . . . . . . . 67ixList of Tables4.1 Summary of simulations carried out with partial surface fields.For simulations 1-10, Lx, Ly, and Lz are 2.67, 10, and 5 nm,respectively, for simulation 11, 2.68, 10, and 7 nm, respec-tively. In all simulations the density is ∼ 0.96 g/cm3. Thenumbers given in brackets in column one are the (y, z) dimen-sions (in nm) of the field bands (Fig. 2.3). The times given inthis table are the lengths of the simulations. For rectanglesthese are the length and width, for triangles the base andheight, and for semicircles the radius. The numbers of cubic,hexagonal, and liquid water molecules in the configuration atthe final timestep are given in the last three columns. . . . . 865.1 Melting temperature estimates for different fields obtainedfrom simulations with 8000 molecules. Estimates of ∆Sfus areobtained from NPT simulations with 1000 molecules. Thelast two columns are the average polarization for liquid andsolid systems. L and S stand for liquid and solid respectively. 1015.2 Summary of ice-like molecules and ice clusters observed indifferent systems; configurations from each simulation weresaved at 10 ps intervals. For some systems, multiple inde-pendent simulations were performed as indicated. At 270 K,freezing occurred within 40 ns for each of the ten simulationsperformed with a field of 1 V/nm; the tabulated data weretaken from pre-freezing configurations. For all systems, thelargest cluster reported is the largest ice cluster (see text fordefinition of cluster) observed that subsequently melted. . . . 105xList of Tables5.3 Summary of ice-like molecules in starting configurations usedto estimate the critical ice nucleus size. The starting config-urations were taken from simulations performed with 32000water molecules interacting with a field of 1 V/nm at 270 K,where freezing eventually occurred. The number of moleculesmeeting the condition for interfacial ice (NI), as identified bythe CHILL algorithm, are recorded in the table, but only cu-bic (NC) and hexagonal (NH) ice-like molecules were includedin determining the size of ice clusters. These four startingconfigurations were used in simulations to estimate the high-est temperature at which the initial ice cluster nucleates ice,both with and without a field. . . . . . . . . . . . . . . . . . 1105.4 Densities in polarized and unpolarized water from simulationswith 8000 molecules at 1 bar. . . . . . . . . . . . . . . . . . . 118xiList of Figures2.1 Schematics of the TIP4P/Ice (left) and six-site (right) watermodels. The water oxygen atoms (O) are red and hydrogenatoms (H) are black, the sites of negative charge on oxygen(M) are cyan, and the lone pair charge sites (LP) are yellow. 172.2 Schematic of the confined water systems. The surfaces aredenoted as two vertical lines, and the water molecules areplaced between them. . . . . . . . . . . . . . . . . . . . . . . 192.3 A sketch of the simulation cell showing the electric field ge-ometries considered. The surfaces are placed at both the leftand right side of the box, and the field is directed along they axis (perpendicular to the page). . . . . . . . . . . . . . . 262.4 A plot showing the trend of the electric field magnitude, wherea = 10 A˚, c = 10 A˚, and f(y, z) = z in Equation (2.15). Thesurfaces are placed at z = 0 A˚ and z = 40 A˚. . . . . . . . . . 273.1 Configurational snapshots of one particular simulation with1200 water molecules. The field (directed upwards along the yaxis) was applied at t = 0 (upper panel), the middle (t = 2.4ns) and bottom (t = 4 ns) panels show the nucleation andgrowth of ice. The oxygen atoms of the water molecules thatexperience the field (z . 10 A˚) are green, those outside thefield region are red, and all hydrogen atoms are black. . . . . 513.2 The dipole order parameter profile along the z axis for thefield-nucleated, frozen sample. The red, green, and blue linesare the x, y, and z components, respectively. Note that they component shows field-induced polarization for z . 10 A˚. . 52xiiList of Figures3.3 Oxygen-oxygen and hydrogen-hydrogen radial distribution func-tions. The black, red, and dark blue lines are the zero field,constant field, and field-nucleated results, respectively. Thelight blue lines represent ferroelectric bulk cubic ice frozenunder a constant field. . . . . . . . . . . . . . . . . . . . . . . 543.4 Configurational snapshots after 18 ns of two selected trial runsthat show ice nucleation and growth. In one system (top),Lz = 60 A˚ (648 particles) and in the other (bottom) Lz = 80A˚ (900 particles). The oxygen atoms of water molecules thatexperience the field are blue, those outside the field regionare red, and all hydrogen atoms are black. The top andbottom projections can be recognized as the (101) and (001)crystallographic planes of cubic ice, respectively. . . . . . . . 573.5 The water-water contribution to the configurational energyfor the six-site model at 270 K. Results for eight trial runs ofthe 432 particle system are shown. The field was appliedat 2 ns. Note that all eight trials show clear nucleationand freezing, as signalled by a rapid drop in the water-waterinteraction energy. . . . . . . . . . . . . . . . . . . . . . . . . 583.6 Oxygen-oxygen radial distribution functions for the six-site(top panel) and TIP4P/Ice (bottom panel) models. Curvesfor different system sizes are shown, and compared with thebulk cubic ice result. . . . . . . . . . . . . . . . . . . . . . . . 603.7 Density (top panel) and dipole order parameter (bottom panel)profiles for two trial runs of the six-site model with 432 particles. 613.8 Configurational snapshots of the two trial runs for whichdensity profiles are shown in Fig. 3.7. Both snapshots areof the configuration at the 32nd nanosecond, and both arefrom the same perspective. The oxygen atoms of the watermolecules that experience the field are blue, those outside thefield region are red, and all hydrogen atoms are black. The icecrystals clearly have different orientations in the simulationcell. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62xiiiList of Figures3.9 Oxygen-oxygen radial distribution functions and mean squaredisplacement (MSD) curves for the six-site model at differenttemperatures. . . . . . . . . . . . . . . . . . . . . . . . . . . . 643.10 Dipole order parameter profiles and oxygen-oxygen radial dis-tribution functions at different values of Emax and c for boththe six-site and TIP4P/Ice models. The values of Emax×10−9V/m and of c are indicated in the figure. . . . . . . . . . . . . 684.1 Configurational snapshots showing ice nucleation and growthfor a 2400 molecule system with a full surface field. Thesnapshots show the (110) plane of cubic ice at 20 ns (stage1), 50 ns (stage 2), 60 ns (stage 3), and 70 ns (stage 4).The water oxygen atoms in the field region are dark blue,those associated with cubic ice or liquid are red, and those inhexagonal ice layers are light blue. All hydrogen atoms areblack. The rectangles outline the central cell. The dashedgreen lines shown in stage 1 indicate growth at (111) planesof cubic ice. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 834.2 Configurational snapshots showing ice nucleation and growthfor a rectangular field band (simulation 2 of Table 4.1). Fromleft to right the snapshots correspond to 15, and 90 ns. Theatoms are coloured as in Fig. 4.1. The dashed green linesin the snapshot at 15 ns indicate ice growth at (111) faces ofcubic ice. Note that the hexagonal ice layers meet the surfaceoutside the field region. . . . . . . . . . . . . . . . . . . . . . 884.3 Configurational snapshot (at 85 ns) showing ice nucleated bya rectangular field band in a larger system (simulation 11 ofTable 4.1). The atoms are coloured as in Fig. 4.1. Note thatin the larger system hexagonal ice layers can form parallel tothe surface. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89xivList of Figures4.4 Configurational snapshots (at 90 ns) of ice nucleated witha semicircular field band (simulation 10 of Table 4.1). The(001) (left panel) and (110) (right panel) planes of cubic iceare shown. The atoms are coloured as in Fig. 4.1. Note thatthe hexagonal ice layers meet the surface outside the fieldregion. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 915.1 Example starting configuration containing ∼ 4000 ice-like molecules.The oxygen atoms of cubic ice-like molecules are shown indark blue, hexagonal ice-like molecules in light blue, andliquid molecules in red. All hydrogen atoms are white. . . . 1025.2 Average number of clusters (per configuration) of a given sizeobserved in NPT (1 bar) simulations with 8000 molecules.Results are shown for: 0 V/nm, 260 K (red); 1 V/nm, 280K (blue); 1 V/nm, 290K (black). Standard errors estimatedfrom ten independent simulations of 40 ns (260 K) or 50 ns(280 and 290 K) are also shown. . . . . . . . . . . . . . . . . 1065.3 Oxygen-oxygen radial distribution functions for ice-like moleculesobtained from simulations (N = 8000, P = 1 bar) for threesystems, as indicated in the legend. . . . . . . . . . . . . . . 1075.4 Example starting configuration with an ice “nucleus” of ∼ 500molecules. Only cubic (blue) and hexagonal (red) ice-likemolecules are shown. . . . . . . . . . . . . . . . . . . . . . . 1085.5 Changes in thermodynamic properties (top) and the dipoleorder parameter (bottom) as functions of field strength. . . . 1125.6 Entropy change as a function of polarization (< mx >) at300 K. The red circles are from Equation (5.13), and the bluetriangles denote the simulation results. The lines are to guidethe eye. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1155.7 Comparison of tetrahedral order parameter distributions atconstant pressure for different temperatures and external fields,as indicated in the legend. . . . . . . . . . . . . . . . . . . . 116xvList of Figures5.8 Total (bottom) and partial (as labeled) pair correlation func-tions at the temperatures and fields indicated in the legend.Multiplication by r2 magnifies the structural features. . . . . 1175.9 Structure factors for the temperatures and electric fields in-dicated in the legend. . . . . . . . . . . . . . . . . . . . . . . 120xviAcknowledgementsForemost, I would like to express my sincere gratitude to my supervisor,Prof. Gren Patey, for all his patient help, enthusiastic discussion, and kindencouragement in my research throughout these years. I also benefit muchmore from his positive and energetic attitude in daily life. Working andtalking with him is very interesting and comfortable, and our conversationis always valuable when I encountered difficulties. I would not have ac-complished my thesis without his guidance, and could not imagine a bettersupervisor for my Ph.D study. I would like to thank all my colleaguesfor their help and advice. When I started my new life and research inCanada, I received lots of support and help from our group members. Ialways appreciate the kindness and warm heart of Dr. Heidrun Spohr,Erin Lindenberg and Dr. Timothe Croteau. I would like to thank Dr.Sarah Overduin and Dr. Debashree Chakraborty for their knowledge andexperiences in simulations and thermodynamics. They were great mentorsin my courses and research, and also good friends for drinking, hiking, anddoing yoga together. Last but no least, I would like to thank my committeemembers and Westgrid facilitators for their help.Finally, I would like to thank all my dear friends for sharing so manyprecious memories with me. I would like to thank my husband Dr. ZhengYang and my parents. I owe a great debt to my family. I would never behere without their love.xviiChapter 1IntroductionWater is one of the most important substances on Earth. It is the vitalsolvent in living cells, the most widely used medium in industry, and an in-tegral part of the hydrologic cycle. Also, water exists as three basic phases onEarth: solid, liquid, and gas, and there are more than 10 phases of ice. Howliquid water turns into ice is crucial in climate change as well as industrialapplications. Why and how ice forms, and what conditions promote icenucleation, are significant questions that remain poorly understood. Theobjective of this thesis is to understand one possible microscopic nucleationmechanism, and the associated thermodynamic properties of the nucleationprocess. As an ideal tool for investigating nucleation, Molecular Dynamics(MD) simulations can provide a direct microscopic picture of ice nucleationand growth.1.1 Nucleation in the Atmosphere and inBiological SystemsIn the hydrologic cycle, water vapour in the air rises mostly by convectioninto the atmosphere, where cooler temperatures cause it to condense intoclouds. Composed of liquid droplets or frozen crystals with water and various11.1. Nucleation in the Atmosphere and in Biological Systemschemicals, clouds have a major effect both short term on the weather, andlong term on the absorption and scattering of solar radiation [1–7].There are two ways for ice to nucleate. Homogeneous ice nucleationproceeds from liquid water in the absence of a foreign substance. In con-trast, heterogeneous nucleation occurs when water is in contact with certainforeign substances, which act as a so-called ice nuclei (IN) [7, 8]. Bothnucleation mechanisms take place in clouds, where ice formation can stronglyimpact the atmospheric conditions and play a major role in precipitation [9].Heterogeneous ice nucleation happens readily in the atmosphere, trig-gered by aerosol particles [3, 10]. The formation of ice clouds is thoughtto be influenced by agents of ice nucleation present in atmospheric aerosols[11]. Atmospheric aerosols are a complex mixture of solid particles andliquid drops, with mineral dust particles being among the most important icenuclei in the atmosphere [12–14]. Large quantities of mineral dust particlesare aerosolized into the atmosphere from various sources, most from aridregions in Africa, the Middle East and Asia [15]. Mica and quartz aremajor contributors to mineral dust, while kaolinite, chlorite, and other clayminerals are found in lower concentrations [16]. These particles provide sitesfor surface reactions and condensed-phase chemistry to take place in theatmosphere, and act as cloud condensation nuclei [6, 13, 17–19]. Kaoliniteparticles are known to be effective ice nuclei and facilitate ice nucleation andcloud formation processes [6, 20, 21]. This could lead to more precipitationand less solar radiation reaching the surface of Earth. Understanding theinteraction of water with mineral dust particles is important for research onthe mechanism of ice nucleation, and indirectly on climate change.21.2. Nucleation TheoryApart from its importance for atmospheric phenomena, the relevanceof ice nucleation is widespread. In astrophysics, abundant subsurface iceis found on many planets in the Solar System [22, 23]. Moreover, icenucleation plays a key role in ecosystems [24] and biological systems [25, 26].Intracellular ice is lethal to cells and in cryobiology it is crucial to control iceformation [27]. Certain epiphytic bacteria capable of living on leaf surfacesare able to catalyse ice formation at temperatures close to 0 ◦C, which resultsin frost injury in plants [28, 29]. Ice accumulation on surfaces is also akey factor contributing to energy efficiency in aircraft de-icing, powerlinemaintenance, and building construction [30, 31]. The applications of icenucleation extend from physiology [32, 33] to food science [34, 35], and fromsnow production [36] to nano-techniques [37].1.2 Nucleation TheoryClassical Nucleation Theory (CNT) [38–41] is a good starting point fromwhich to consider the ice nucleation process. CNT uses macroscopic ther-modynamic quantities to describe the emergence of a new phase. It can beused to describe both homogeneous and heterogeneous nucleation processes.Homogeneous nucleation refers to the creation of a new phase in the absenceof any foreign substance, while heterogeneous nucleation is facilitated by theaction of foreign agents or surfaces.Near the equilibrium transition point, the original phase remains metastable,and the formation of a new phase depends on the height of the Gibbs freeenergy barrier. At temperatures below 0 ◦C, supercooled liquid water is31.2. Nucleation Theorymetastable with respect to natural crystalline ice, and can form new solidclusters. The assumption is made that the solid cluster is spherical and hasa radius r, thus the free energy difference between the two phases can bewritten as [1–3, 42]∆Gp =4pir33v (µ2 − µ1), (1.1)where the subscript p indicates the free energy change due to the phasechange, v is the molecular volume of ice, and µi is the chemical potentialof phase i. This term is always negative and favourable, but creating aninterface between the two phases is associated with an increase in free energy∆Gs = 4pir2σ, (1.2)where the subscript s indicates the free energy change due to surface tension,σ is the surface tension between the new and parent phases. The nucleationprocess is controlled by the overall Gibbs free energy change, which is thesum of the volume free energy change of Equation (1.1) plus the surface freeenergy of Equation (1.2)∆Ghomo =4pir33v (µ2 − µ1) + 4pir2σ, (1.3)where the subscript homo indicates the free energy change for a homoge-neous nucleation process.As the cluster size becomes larger, the overall free energy change willincrease and reach a maximum, which is called the free energy barrier tonucleation. At larger sizes, the free energy will drop and lead to ice growth.41.2. Nucleation TheoryThe position of the free energy maximum determines the critical radius r∗,which can be obtained by taking a derivative [1, 42] of Equation (1.3)∂∆Ghomo∂rr=r∗= 0. (1.4)The critical radius r∗ and the critical free energy ∆G∗homo arer∗ = −2σv(µ2 − µ1), (1.5)∆G∗homo =16piσ3v23(µ2 − µ1). (1.6)When new born clusters are smaller than the critical size r∗, a moleculeleaving the clusters is preferable to a molecule being added. Therefore, icenucleation will not be favoured and clusters may decompose back to liquid.However, for clusters larger than r∗, growth will be spontaneous by conden-sation of further molecules, and the free energy decreases continuously.Usually, the critical conditions required for homogeneous nucleation canbe quite difficult to achieve in real systems. However, the concepts andprinciples of homogeneous nucleation are simple and helpful for analysingother types of nucleation. The equation of critical free energy [Equation(1.6)] highlights the strong dependence on interfacial tension, and it is easyto speculate that a suitable foreign substances in contact with liquid watercan reduce the free energy barrier and catalyze ice nucleation.Let us assume that heterogeneous nucleation is taking place on an insol-uble foreign flat surface, and the newly formed ice cluster is hemisphericalon the planar surface. In this case, the classical expression for the critical51.3. Field Induced Ice Nucleationfree energy ∆G∗het of cluster formation on the surface will be [1, 2]∆G∗het = ∆G∗homof(θ), (1.7)where θ is the contact angle between the spherical cluster and the flatsubstrate, and f(θ) is a geometric factor with values between 0 and 1. Bythe presence of a surface, the height of the free energy barrier relative tohomogeneous nucleation is reduced by the factor f(θ). Although f(θ) canbe established empirically, its physical significance is unclear since the icenucleus may not adopt a hemispherical form [1, 3]. In some cases f(θ) maycorrespond to adaptations in a single molecular layer at the interface betweenthe surface and the new cluster, whereas in other cases the reduction of thefree energy barrier by a substrate may be due to density fluctuations causedby the impact of the substrate. Since pure bulk water is rare in daily life,and heterogeneous nucleation is effective at warmer temperatures [43], theformation of water droplets and ice in nature is usually by heterogeneousnucleation.Classical nucleation theory assumes a spherical nucleus, and oversimpli-fies the properties of clusters and surfaces in ice nucleation, therefore onemust be aware of its limitations [42].1.3 Field Induced Ice NucleationWhat promotes ice nucleation has been an interesting topic for researchersfor years, and both experimental and theoretical methods are employed to61.3. Field Induced Ice Nucleationreveal the nature of the ice nucleation process. Although scientists havenoticed connections between electric fields and ice formation, electrofreezingof bulk water has not been firmly established experimentally [44, 45], butelectrically induced or assisted freezing in particular environments has beenreported. Bulk water undergoes dielectric breakdown under uniform fieldsat magnitudes of 6.5 - 7.0 ×107 V/m [46]. As pointed out in Ref. [44], theelectric fields achievable on macroscopic length scales are orders of magni-tude smaller than the field strength (∼ 109 V/m) required to electrofreezemodel water in computer simulations [47]. For computational simulations,the high magnitude of the field needed to nucleate ice is perhaps due tothe time limitation of simulations, and a lower field might work on muchlonger time scales. However, fields of comparable strength to theoreticalprediction (∼ 109 V/m) do exist on microscopic length scales, leaving openthe possibility of field-induced ice nucleation.There have been a number of experimental reports of ice nucleation inthe presence of strong electric fields. The results are in general positiveconcerning the influence of field on ice nucleation, but remain conflicting.In the 1960s, it was reported that highly charged silver iodide particlesprohibited ice nucleation on the surface [48]. Later Dawson and Cardell[49] observed that an electric field did not increase the freezing possibilityof water at temperatures between −8 ◦ and −15 ◦C. Doolittle and Vali[50] applied an electric field of up to 6 × 105 V/m over sets of super-cooled drops supported on a flat solid surface, and found that there wasno significant difference in the heterogeneous nucleation rate due to anelectric field. On the other hand, Gabarashvili and Gliki [51] found that71.3. Field Induced Ice Nucleationsupercooled water drops containing impurities with charged surfaces werenucleated to ice at significantly warmer temperatures when the particlescarried a net negative charge than when they were neutral or carried anet positive charge. It was also found by Abbas and Latham [52] thatthe probability of supercooled drops to freeze was greater if the surface ofthe drop was disrupted by electrical forces to produce a filament of liquidfrom the localized area of rupture. Later, Pruppacher [53] found that thefreezing temperature of supercooled water was raised considerably when incontact with predominantly negatively charged amorphous sulfur particles,and concluded that electrostatically charged surfaces and local electric fieldssignificantly enhance ice nucleation. The decadal data from 1953 to 1985on high-energy cosmic ray fluxes indicates that the atmospheric ionizationaffects ice nucleation in high-level clouds, and the electrofreezing increasesthe flux of ice crystals [54].The work of Garvish et al. [55] represents a new milestone in the researchof ice nucleation under local electric fields. There authors reported that polaramino acid crystals induced a freezing point higher by 4 - 5 ◦C than non-polar crystals. These results were interpreted in terms of an electric fieldmechanism that helps align the water molecules into a proton ordered ice nu-cleus. Braslavsky and Lipson [56] have reported that a high-voltage electricpulse applied to bulk supercooled water induces ice nucleation. Also, Ehreand coworkers [57] have reported that water freezes differently on positivelyand negatively charged surfaces. A surface with positive charge promotedice nucleation, whereas the same surface when negatively charged reducedthe freezing temperature. This fact may be related to the orientation of81.3. Field Induced Ice Nucleationwater molecules with respect to the surface, but the details of how thiseffect occurs are not yet understood.The situation for computer simulations of water subjected to electricfields is clearer. Svishchev and Kusalik [47, 58] were the first to show, em-ploying MD simulations, that bulk water freezes readily to form ferroelectriccubic ice Ic when polarized by a strong uniform external field (5×109 V/m).The models used were the TIP4P [59] and the SPC/E [60] water models, andthe former appeared to be more susceptible to ice nucleation than the latter.Other authors have demonstrated electrofreezing of water confined betweenplates. Xia and Berkowitz [61] found that when the charge density on Ptsurfaces reaches a threshold value, water between these surfaces can freezeinto cubic ice. Zangi and Mark [62] reported that water confined betweentwo plates freezes into ice when an electric field is applied parallel to theplates. While these studies suggest that electric fields might be importantfor ice nucleation, they do not directly address the role of electric field inthe electrofreezing process.Although the field strength necessary to cause freezing in simulations(∼ 109 Vm−1) [47, 62] is much larger than fields that can be establishedin macroscopic samples [44], they may exist on smaller length scales. Forexample, a large electric field can occur along the axis of a polar amino acidcrystal [55]. Also, for trench-like structures in model kaolinite, there is anelectric field due to surface charge arrangements, which can strongly polarizeand order water [63]. Local electric fields can exist in nature, and couldpossibly be a key factor in ice nucleation. There are, however, importantquestions to be solved, such as how an ice nucleation process is triggered by91.4. Thesis Motivation and Objectivesthe field, and the initial structure and growth of ice nuclei.1.4 Thesis Motivation and ObjectivesDespite the importance of ice nucleation discussed above, the simple ques-tion: What makes a good IN? has proved difficult to answer. One often citedexplanation invokes a good crystallographic match between an IN surfaceand hexagonal ice [7]. Kaolinite is sometimes cited as an example, becausekaolinite particles are good IN, and one of its basal planes (the Al-surface)has a potential crystallographic match with hexagonal ice. However, thisidea is not supported by recent model calculations [64–67]. A classicalmodel [64] shows that water adsorbed on this surface does not take on anice-like structure, and coverage beyond a bilayer is not found. Therefore, thecrystallographic match argument does not appear to explain kaolinite’s icenucleating abilities. A similar conclusion was reached for model hexagonalsurfaces specifically selected to provide a good match with hexagonal iceIh [65]. Moreover, an extensive density-functional theory (DFT) study [67]shows that water can wet kaolinite’s basal surface, but multilayer ice growthis not favoured. Hu and Michaelides [66] suggest that the stable overlayerof water on a kaolinite surface is not a consequence of the lattice matchwith Ih, but rather because the substrate is amphoteric with the ability toform both acceptor and donor hydrogen bonds. A recent MD simulationof ice nucleation near a kaolinite slab [21] shows that on the basal face ofkaolinite, growth along the prism face of ice is observed, which contradictsthe expectation of the lattice match between ice and the kaolinite surface.101.4. Thesis Motivation and ObjectivesExperiments [68] tend to support the idea that ice nucleation occurs at“active sites”, but the nature of such sites remains an open question.Studies on kaolinite reveal that the effectiveness of kaolinite for ice nu-cleation is likely more interesting and complex than a simple lattice match.A quantum mechanical calculation [69] has demonstrated that the kaolinitebasal surface is polar. Croteau et al. have shown [63] that for particular,trench-like structures in kaolinite, the field due to surface charge arrange-ments can strongly polarize and order water over distances correspondingto a few molecular layers. An interesting MD study of homogeneous icenucleation reported by Vrbka and Jungwirth [70] suggested that the electricfield associated with water ordering at the surface might partially explain theobservation that ice preferred to nucleate in subsurface layers. Note that thiswork could not be repeated for large systems perhaps due to computationallimitations [70]. Inspired by these works, we consider a model of confinedwater experiencing a surface electric field, which polarizes a thin layer ofwater molecules near the surface. Note that in this model, the surface fieldonly exists near the surface, and decays to insignificance very rapidly withdistance from the surface.As noted above, it has long been suggested [53] that electric fields couldpossibly promote freezing, and electrically induced freezing has been con-firmed by simulations [47, 62]. However, computational simulation studieshave usually focused on electrofreezing, which refers to the process where auniform electric field is applied to the whole system, and can change thermo-dynamic properties, such as the melting point itself. Although electrofreez-ing has been established, there are no previous simulation studies showing111.4. Thesis Motivation and Objectivesice nucleation induced by a partial field, where the field only impacts a smallpart of the system. In this thesis, we demonstrate a different mechanismfor heterogeneous ice nucleation when water interacts with a surface field.Our simulation results confirm the speculation that a microscopic surfacefield can have a a positive influence on ice nucleation. We show that waterfreezes into a thin polarized layer of ice with the help of the surface field,and then grows further beyond the range of the field to form bulk ice.The overall goal of this thesis is to use computer simulations to inves-tigate the microscopic mechanism of ice nucleation induced by an electricfield, and to understand the patten of ice growth. Firstly, a simple setupis considered, where water molecules located very near a surface (within∼ 10 A˚) experience an electric field parallel to that surface. The surfacesconsidered are more or less “neutral”, and do not induce any significantorder in the water. In supercooled water, we show how water moleculesinteract with the surface field, and and how ice nucleation can occur. Thestructure of the system, evolving in time, is analyzed. To ensure thatour observations are not an artifact of a particular water model or choiceof system parameters, simulations are carried out for a range of systemsizes and dimensions. Both the TIP4P/Ice [71] and the six-site [72] watermodels are considered and the results are compared. We also examine theimportant influences of temperature and electric field parameters (strengthand extent) on ice nucleation and growth. Density and polarization profiles,as well as mean square displacement profiles are analyzed in the process ofice nucleation and growth. This study provides insight into some possiblecharacteristics of good IN.121.4. Thesis Motivation and ObjectivesWe further explore ice nucleation by modifying the surface field. Weconsider surface fields which act only over surface “bands” of different sizeand geometry. The same two plates are used, and the electric field is appliedover a particular geometry near one plate, and decays outside the designatedarea. By varying the field geometry, we learn more about the mechanism ofice nucleation. In systems with an electric field band, ice freezes inside thefield region, and eventually grows further beyond the range of the field. Thesize of the ice nuclei induced by the electric field is crucial for ice nucleation.We also look into the pattern of ice growth on different crystallographicplanes. The water-ice structure is followed in time, and the results show apatten of layer-by-layer ice growth on the (111) crystal plane. Comparisonsof ice growth are made between infinite surface fields and surface field bands.Structures of hexagonal ice are very rare in the former system, which is dueto the infinite layer of ice Ic on the surface.Both hexagonal ice and cubic ice can be constructed by stacking bilayersof water molecules formed of hexagonal rings. But how the bilayers arestacked is different in the two ice phases. In the basil planes of ice Ih, twohexagonal layers are stacked following the ABABAB... sequence, while the(111) planes of ice Ic consist of a stack of the same hexagonal double layersbut following the ABCABC... sequence. Then the basal planes of ice Ih andthe (111) planes of ice Ic are identical and can be perfectly joined along theirhexagonal double layers with very small interfacial energy [73, 74]. Thus itis possible to form either structure when ice is growing on the hexagonalplanes.Ice structures are analyzed using the CHILL algorithm originally devel-131.4. Thesis Motivation and Objectivesoped by Moore et al. [75], and home-written code. As hexagonal and cubicice are both made of similar layers of water molecules and differ in theirstacking pattern, a mixture of ice Ih and Ic is possible when ice grows layerby layer. The question as to whether or not the formation of a subsequentlayer on a current (111) crystal layer is random or biased is interesting,and decides the ratio of these two ice structures. The ratio of the two icestructures is reviewed and compared with other work.In addition to studying the phenomenon of heterogeneous ice nucleationinduced by varying fields, we perform simulations on homogeneous ice nu-cleation in the presence of a uniform electric field. The isothermal isobaric(NPT ) ensemble is employed. The goal here is to find the intrinsic reasonwhy electric fields can catalyze ice nucleation, and how water structure isinfluenced by the field. We investigate the structural and thermodynamicchanges that promote freezing when model liquid water interacts with auniform, external field. It is found that the melting temperature increasessignificantly when electric fields of ∼ 1 × 109 V/m are applied; however,liquid water maintains its spatial structure and tetrahedral order, both aboveand below the melting temperature. Although the melting temperatureincreases, water interacting with a field remains liquid until it is ∼ 40 Kbelow the field-dependent melting temperature. The degree of supercoolingdetermines the size of the critical nucleus, regardless of field. Our resultssuggest that electric fields produced by surface charges promote freezing bylocally “supercooling” water molecules, allowing ice clusters to form andgrow.141.5. Thesis Overview1.5 Thesis OverviewThis thesis is presented in six chapters. The current chapter reviewedbackground information relevant to ice nucleation applications, classicalnucleation theory, as well as the current knowledge of electrofreezing andice nucleation. Chapter 2 describes in detail the water models used, simu-lation methods employed, and the methods of analysis. Chapter 3 discussesheterogeneous ice nucleation induced by varying surface fields in confinedsystems. Chapter 4 deals with ice nucleation induced by field “bands” ofdifferent geometries. Chapter 5 considers the mechanism of electrofreezing,and finally, a summary of the results and future directions are presented inChapter 6.15Chapter 2Models and SimulationMethodsIn the process of a molecular dynamics (MD) simulation, potential energyand kinetic energy are calculated by knowing the configuration at each timestep, and the system propagates by solving the equations of motion. Thischapter describes models of the system with water and surfaces, and themethods used in MD simulations. The code that I wrote is compared withopen source packages, and the data analysis is discussed.2.1 Molecular ModelsThis thesis uses rigid water models to investigate field induced ice nucleationprocesses. We consider the six-site [72] and the TIP4P/Ice [71] water models.As shown in Fig. 2.1, both water models have an oxygen atom (labeled O) atthe center, hydrogen atoms (labeled H) with positive charges on both sides,and also a “dummy” atom (labeled M) with negative charge sitting near theoxygen, on the bisector of the HOH angle. Differing from the TIP4P/Icemodel, the six-site model places negative charges on two additional dummyatoms (labeled LP) representing the “lone pairs” electrons of the oxygen162.1. Molecular Modelsatom. Parameters of both water models are given in Table 2.1.Figure 2.1: Schematics of the TIP4P/Ice (left) and six-site (right) watermodels. The water oxygen atoms (O) are red and hydrogen atoms (H) areblack, the sites of negative charge on oxygen (M) are cyan, and the lone paircharge sites (LP) are yellow.For confined water systems, stationary surfaces (labeled S) are placedas shown in Fig 2.2. The surfaces are considered as smooth and solid, andinteract with the oxygen atom of each molecule through the potentialu(z) = 3√3ε2[(σz)9−(σz)3], (2.1)where z is the prependicular distance of the oxygen atom from the surface,and ε and σ (given in Table 2.1) are the energy and distance parametersof the surface. This potential has been previously used to model smoothsurfaces that interact weakly with water and impose little order in the liquid172.1. Molecular ModelsParameter Six-site TIP4P/Ice Surfaceσ (A˚) 3.115OO 3.1668OO 2.473OS0.673HHε (KJ/mol) 0.715OO 0.8822OO 1.932OS0.115HHr1(A˚) 0.980 0.9572r2 (A˚)0.230OM 0.1577OM0.8892OLq+ (e) +0.477 +0.5897q− (e) −0.866M −1.1794M−0.044Lθ(deg) 108.00HOH 104.52HOHφ(deg) 54.00MOH 52.26MOH111.00LOLTable 2.1: Water model and surface parameters used in the simulations.Note that ε and σ are the energy and distance parameters for the interactioncalculation, r1 is the oxygen-hydrogen bond distance, r2 is the distancebetween the oxygen and the negative charges (OM or OL), q is the chargeon each site, and θ and φ are angles associated with different sites in watermolecules. The subscripts HOH, MOH and LOL indicate the angles betweentwo hydrogen atoms, between the negative charge M and the hydrogen atom,and between the two lone pair charge sites. Note that L is short for lonepair (LP).[76].For all systems considered, periodic boundary conditions are implementedin order to eliminate surface effects. The original simulation box is replicatedthroughout space to form an infinite lattice. The periodicity retains a central182.1. Molecular ModelsFigure 2.2: Schematic of the confined water systems. The surfaces aredenoted as two vertical lines, and the water molecules are placed betweenthem.simulation cell with a constant number of molecules, and also allows watermolecules to interact with each other and with periodic images. Bulksystems are periodic in x, y and z directions, and for confined systems,we adopt a simulation cell consisting of two parallel walls and two emptyregions outside of the walls, which is replicated in the x, y and z directions.The empty area should be large enough to get extremely small interactionsbetween water confined in the original box and the periodic boxes in the zdirection.Water molecules interact with each other through Lennard-Jones (LJ)potentials between atoms and electrostatic (Coulombic) interactions be-tween charges [77]. In the TIP4P/Ice water model, the LJ potential betweendifferent molecules is only located on the oxygen atoms, while in the six-sitemodel, both oxygen and hydrogen atoms interact with each other throughLJ potentials. In both models, charged sites account for the electrostaticpotential. The configurational energy of water is calculated by summing up192.2. Interaction Potentialsthe Lennard-Jones and Coulombic interactions over all sites in all watermolecules. There is an additional water-wall potential for the confinedsystems when surfaces are present. Details of potential calculations arediscussed in the next section.2.2 Interaction PotentialsAs noted above, the intermolecular interactions are modelled through Lennard-Jones potentials between pairs of atoms, and Coulombic potentials actingbetween pairs of charges such thatu(rij) = 4ε[(σijrij)12−(σijrij)6]+ uij,Coul, (2.2)where rij is the distance between two atoms, εij and σij are the energyand distance parameters. Note that the second term in Equation (2.2) onlyaccounts for sites with charges. The LJ potential is truncated spherically ata reasonable cutoff distance (∼ 10 A˚) to improve the efficiency of calculation,which may vary a little in different systems. For interactions betweendifferent types of atoms, LJ parameters are calculated according to theLorentz-Berthelot mixing rules [78, 79]σij = (σi + σj)/2, (2.3)εij =√εiεj , (2.4)202.2. Interaction Potentialswhere σi and εi are listed for the different models in Table 2.1. The totalLJ contribution to the configurational energy is a summation of pairwise LJinteractions over all pairs of atoms.In our model, there are Coulombic interactions between electrically chargedatoms or sites in different molecules. The Coulombic potential is given byuij,Coul =14piε0qiqjrij(2.5)where ε0 is the permittivity of free space, qi and qj are a pair of charges iand j, and rij is the distance between the charges.The Coulombic potential is a long-range interaction. The long-rangeinteraction falls off slowly as 1/rij , and its significant range is usually greaterthan half the simulation box length. The charges on each molecule will notonly interact with molecules in the central simulation cell, but also with allimages in the periodic cells. The total Coulombic potential energy of thesystem can be expressed asUCoul =18piε0∑nN∑i=1N∑j=1′ qiqj|rij +nL|, (2.6)where L is the original box length, N is the total number of charges, rijis a vector joining site i and j, and n is an integer vector denoting eitherthe original simulation box or periodic image. The prime indicates thatwe omit i = j for n = 0. Further periodic images are constructed and asphere of simulation boxes is generated, such that interactions with ions inthe original box can be explicitly calculated [77].212.2. Interaction PotentialsIn practice, Equation (2.6) is slowly converging and conditionally con-vergent. The Ewald method can be used to tackle the long-range forceswith efficiency, and was developed to sum the interactions between chargesand all of their periodic images [77, 80, 81]. The charge density distributionof each point charge in the system is split into two terms by adding andsubtracting a Gaussian distribution. It is assumed that each point chargeis surrounded by a Gaussian charge distribution of different sign, with thefunctional formρi(r) = qiκ3exp(−κ2r2)/pi3/2, (2.7)where the parameter κ determines the width of the Gaussian distribution,and r is the position relative to the center of the distribution. This extradistribution serves to screen the interaction between point charges, suchthat the screened interactions are now short-ranged and will converge morerapidly. A cancelling Gaussian distribution with the opposite sigh of Equa-tion (2.7) is considered. This cancelling distribution recovers the overallpotential of the original set of charges. By taking a Fourier transform, thecancelling distribution can be summed in reciprocal space.In other words, Equation (2.6) is split into two terms, which are calledreal space (Ur) and reciprocal space (Uk) terms. However, all charges aretreated equally in the Ewald summation technique, and the interactionsbetween charges from the same molecule are included. In order to eliminatethe unwanted interactions between charges within the same molecule, twoadditional terms must be subtracted. Thus, the total Coulombic potential222.2. Interaction Potentialsenergy isUCoul = Ur + Uk − Us − Ud, (2.8)where Ur, Uk, Us, Ud are the real space term, the reciprocal space term, theGaussian term, and the self-energy correction, respectively. Each term isdescribed in detail as follows.In the Ewald method, the real space energy calculation is truncated bythe erfc(x) function asUr =18piε0N∑iN∑j∑n′ qiqj erfc(κ|rij + n|)|rij + n|, (2.9)where erfc(x) is the complementary error function. κ is usually chosen suchthat only (n = 0) contributes to UrThe summation calculated in reciprocal space isUk =12ε0L3∑k 6=0N∑iN∑jqiqjk2 exp(−k2/4κ2)cos(k · rij), (2.10)where k = 2pin/L2 is a reciprocal space vector. For systems with hundreds ofwater molecules, we typically choose κ = 6/L for the real space calculation,and 1200 k-vectors in reciprocal space. The real space contribution willbe negligible beyond some cutoff, and is confined to the central simulationbox (n = 0), and the reciprocal term can be easily calculated with 1200reciprocal lattice vectors.In addition, a self-energy term, Us, for interactions of the cancellingGaussian distribution with itself should be subtracted. For molecular sys-232.3. Models of Electric Fieldstems, another self-energy term, Ud, accounting for intra-molecular inter-actions should also be taken out of the total potential energy [77]. TheGaussian term Us isUs =κ4pi3/2ε0N∑iq2i , (2.11)where the summation is over all charges qi. The self-energy term for amolecular system isUd =14piε0∑i(∑aqa∑bqberf(κ|d|)|d|), (2.12)where a and b are different charges in the same molecule i and d is the vectorbetween a and b [77]. Note that Equation (2.8) allows for the calculationof the Coulombic potential energy of a simulation cell surrounded by aconductor with a dielectric constant ε′ = ∞.Forces due to Coulombic interactions are calculated for each term inEquation (2.8).2.3 Models of Electric FieldsA polar liquid, such as water, will undergo reorientation when exposed toan external electric field. The field interacts with the molecule by exertinga force on each charge. A uniform field will not influence the translationalmotion of the molecule, but rotate the dipole to the field direction. Thewater-field interaction energy is calculated through a summation over all242.3. Models of Electric Fieldscharges in the systemUef = −N∑i=1qi(ri ·E), (2.13)where qi and ri is the magnitude and the position vector of the ith charge,N is the total number of charges in the system, and E is the electric fieldvector. The force on the i-th charge qi due to the external field is calculatedbyFief = qi · E. (2.14)In this thesis, confined water systems are investigated, where moleculesnear one surface experience an external electric field parallel to that surface.The simple model we first consider is a decaying field which extends somedistance from the surface, and dies off very quickly. The half value of thefield is marked as a yellow line in Fig. 2.3. Also, we consider field bands(infinite in x, and finite in y and z) with cross sectional areas that arerectangular, triangular, and semicircular in shape, as shown in Fig. 2.3.This electric field retains its maximum magnitude inside the field bands,and then decreases rapidly to zero.The functional form selected for the field model is rather arbitrary. Thefield E(z) is given by the functionE(z) = Emax( 11 + ea(f(y,z)−c))yˆ , (2.15)where Emax is the maximum magnitude of the field, f(y, z) is a functionof the molecular position relative to the surface, yˆ is a unit vector alongthe y axis, and a and c are parameters determining the extent and decay252.3. Models of Electric FieldsZXFigure 2.3: A sketch of the simulation cell showing the electric fieldgeometries considered. The surfaces are placed at both the left and rightside of the box, and the field is directed along the y axis (perpendicular tothe page).262.3. Models of Electric Fieldsrate of the field as a function of z. In our model, the field is infinite inx direction, and determined by the y, z coordinates of water. Appropriatechoices of f(y, z) control the shape and extent of the field bands. Note thatthe exponential form ensures a rapid decay as indicated in Fig. 2.3.−1 0 1 2 3 4 5 6 0  5  10  15  20  25  30  35  40Field magnitude ( × 109  V/m)z (Å)Figure 2.4: A plot showing the trend of the electric field magnitude, wherea = 10 A˚, c = 10 A˚, and f(y, z) = z in Equation (2.15). The surfaces areplaced at z = 0 A˚ and z = 40 A˚.For a surface field decaying only in the z direction, we choose a = 10 A˚,c = 10 A˚, and f(y, z) = z, which makes the field fall very rapidly and becomeinsignificant beyond ∼ 10 A˚. The magnitude of the field is plotted as the redline in Fig. 2.4. We emphasize that Equation (2.15) is simply a convenientfunctional form for a field that decays with distance from a surface, it is272.3. Models of Electric Fieldsnot meant to accurately represent any particular physical situation. Inphysical systems, the field near a surface might decay for various reasons,such as screening by mobile charges, or simply from the particular geom-etry and charge distribution associated with electrically neutral surfaces,as in kaolinite trenches [63]. In this thesis, Emax is usually set around5 × 109 Vm−1, which is similar to the value used in earlier constant fieldsimulations [47, 58, 62]. This field is too large to be realistic for macroscopicsystems, but is not inconsistent with fields that can exist on microscopiclength scales (note that the field 6 A˚ from an elementary point charge is4× 109 Vm−1).Generally, the field is kept constant for all point charges located in thesame water molecule. Thus the field exerts a torque but no translationalforce on a water molecule. In real physical situations, we would expect aspatially varying field to exert some force that could influence the densityand packing of water molecules near a surface, and this could vary widelydepending on the local features of a rough surface. However, since the detailsof the surface variation implied by Equation (2.15) are arbitrary, here weinclude only the polarizing aspects of the field, which we would expect to bealways present and less sensitive to details of the field gradient. Note thereare some exceptions with the field model with open source packages, whichwill be discussed later.282.4. System Evolution2.4 System EvolutionIn molecular simulations, the system configurations are evolved in time bycalculating the forces on each molecule and then integrating the equationsof motion. Evaluation of the forces is the most time consuming step in amolecular dynamics simulation. Forces on each object are calculated fromthe interaction potential u(rij) described in Section 2.2,F = −∇u(rij), (2.16)where the ∇ is the gradient operator. For each simulation loop, forces oneach molecule are updated to generate the time trajectory of the systemevolution. In the following sections, how we deal with translational androtational motion will be described separately in detail.2.4.1 The Gear Predictor-Corrector Method forTranslational MotionFor our rigid models, the molecular motion is divided into two parts; thetranslational motion of the center of mass and the rotational motion aboutthe center of mass. Both depend on forces exerted on each molecule.The force experienced by a molecule is the sum of the forces on eachsite within the molecule. From Newton’s equation of motion, the actualacceleration aC of the center of mass is given byaC = Fm, (2.17)292.4. System Evolutionwhere F and m are the total force on and the mass of the particle, respec-tively, and the superscript C refer to the actual value of the acceleration.For translation, the molecule is treated as a whole, and the coordinatesand motional properties are for the center of mass. To begin a simulationloop, the molecular position vector r is estimated from the velocity, accel-eration, and higher derivatives [84].rP (t+∆t) = r(t)+ drdt∆t+12!d2rdt2 ∆t2+ 13!d3rdt3 ∆t3+ 14!d4rdt4 ∆t4+. . . , (2.18)where the the superscript P refer to the predicted value of the positionvectors. Similarly, the particle velocities v and accelerations a can bepredicted asvP (t+ ∆t) = v(t) + d2rdt2 ∆t+12!d3rdt3 ∆t2 + . . .aP (t + ∆t) = a(t) + d3rdt3 ∆t+ . . . ,(2.19)as well as higher derivatives. Predicted values from Equation (2.18) andEquation (2.19) can be used for energy and force calculations, and the actualacceleration aC can be obtained by Equation (2.17). Due to truncationerrors, the predicted values of position vectors and their higher derivativesmust be corrected to yield the right trajectory. The acceleration differencebetween the predicted value aP and the actual aC is∆a(t+ ∆t) = aC(t + ∆t)− aP (t+ ∆t), (2.20)302.4. System Evolutionand this is used to correct all predicted values byrC(t+ ∆t) = rP (t+ ∆t)− k0∆a(t+ ∆t)vC(t+ ∆t) = vP (t+ ∆t)− k1∆a(t+ ∆t),(2.21)where k0 and k1 are Gear corrector coefficients associated with position,velocity, and other derivatives, as listed in Ref. [77]. We calculated upto the 4th derivatives of coordinates in our program. Equation (2.17)is a second order equation, so the Gear corrector coefficients are 19/120,3/4, 1, 1/2, 1/12 for coordinates, velocities, accelerations, and 3rd and4th derivatives of coordinates, respectively. Coordinates and their higherderivatives are corrected as in Equation (2.21) and the proper configurationof the system is obtained. The Gear predictor-corrector method can beviewed as a three-stage procedure. Initially the coordinates and other timederivatives are predicted from the Taylor expansion with a truncation. Thenactual accelerations are determined by new forces for the new configuration.Subsequently, all predicted values are corrected by the acceleration difference∆a(t+ ∆t), and saved for calculation in the next step [84].The Gear predictor-corrector method is not only suitable for transla-tional motion of molecules, but is also a good tool for rotational motion.How the rigid water models are treated for rotation is discussed in the nextsection.312.4. System Evolution2.4.2 Quaternions and Rotational MotionIn simulations, a molecule is moved and rotates around its own center ofmass. In other words, its center of mass is tracked in a global space-fixed coordinate system, while in each molecule, positions of off-center sitesare defined with respect to its own body-fixed coordinate system. Forsmall, rigid, non-linear models, quaternions are introduced to represent theorientation of the body-fixed frame with respect to the global space of thesimulation box [77, 82, 83]. A quaternion parameter Q of a molecule is aset of four scalar quantities [77, 84],Q = (qw, qx, qy, qz). (2.22)The quaternions are not independent but satisfy the relationq2w + q2x + q2y + q2z = 1. (2.23)Quaternions, associated with each molecule, are randomly assigned at thebeginning of a simulation, and the values of (qw, qx, qy, qz) should be rescaledaccording to the magnitude√q2w + q2x + q2y + q2z to satisfy Equation (2.23).For each molecule, it is requisite to know its position and orientation inorder to evaluate all interactions with other molecules. To make the forcecalculations easier, global space-fixed positions of all molecular sites arecrucial. With the help of a rotation matrix A, the space-fixed coordinatesr and body-fixed coordinates d are connected. The space-fixed positions of322.4. System Evolutioneach site a in a molecule i can be obtained by [77, 82, 84, 85]ria = ri +A−1dia, (2.24)where the rotation matrix A is defined as [77, 84]A =q2w + q2x − q2y − q2z 2(qxqy + qwqz) 2(qxqz − qwqy)2(qxqy − qwqz) q2w − q2x + q2y − q2z 2(qxqy + qwqx)2(qxqz + qwqy) 2(qxqy + qwqx) q2w − q2x − q2y + q2z. (2.25)The quaternions also satisfy the equations of motion [77, 82]q˙wq˙xq˙yq˙z= 12qw −qx −qy −qzqx qw −qz qyqy qz qw −qxqz −qy qx qw0ωbxωbyωbz, (2.26)where ωb represents the angular velocity for each molecule in the body-fixedframe.The rotational motion is governed by the torque τ . The torque for eachmolecule can be calculated as the cross product of the forces on each siteand the distance to the center [77, 82],τi =∑adia × fia, (2.27)where dia is the position vector of atom a relative to the center of themolecule i. With the help of the rotation matrix, the torque in the body-332.4. System Evolutionfixed system isτ bi = Aτi = A∑adia × fia, (2.28)and is used to calculate the first derivatives of the rotational velocitiesω˙bx = τbxIxx +(Iyy−IzzIxx)ωbyωbz,ω˙by =τbyIyy +(Izz−IxxIyy)ωbzωbx,ω˙bz = τbzIzz +(Ixx−IyyIzz)ωbxωby,(2.29)where Ixx, Iyy and Izz are the three principal moments of inertia for eachmodel. Note that Ixx, Iyy and Izz are calculated for each water model, anddiffer from the values of real water because of the positions of the atoms.As for the translational motion in Section. 2.4.1, both the quater-nions and angular velocities can be handled by the Gear predictor-correctormethod. Similarly, up to 4th order derivatives are tracked, but the quater-nions and angular velocities follow first-order differential equations, so theGear corrector coefficients [77] for quaternions or angular velocities and theirhigher derivatives are 251/720, 11/12, 1, 1/3, 1/24, respectively.At the beginning of a step in a simulation, the quaternions and angularvelocities are predicted, as are their higher derivatives. Predicted quater-nions associated with positions of the centres of mass are used to obtainthe space-fixed coordinates, and then to calculate interactions and forcesacting on the sites. From the force calculations, the torques are obtained byEquation (2.27), and then the first derivatives of angular velocities can beevaluated according to Equation (2.29). The difference between predictedvalues and values from Equation (2.29) are used to correct angular velocities342.4. System Evolutionand their higher derivatives. Quaternions are treated in the same way asangular velocities. Note that the magnitude of a quaternion must be checkedand rescaled to 1 at every time step.For both translation and rotation, the Gear predictor-corrector method isfast and tolerant of relatively large time steps, for example, 2 femtoseconds.It is one of the most widely used predictor-corrector algorithms for moleculardynamics. The prediction part is applied at the beginning of every time step,and the correction part is applied after the force evaluation. This methodis usually implemented to generate long time simulation trajectories.2.4.3 Simulation ConstraintsMost simulations carried out in this work are in the canonical ensemble,which has a constant number of particles (N), volume (V ) and temperature(T ). A variety of thermostating methods are available to add or removeenergy from the MD system, approximating a constant temperature.The simplest method to control temperature involves velocity scaling[84] asvnew = vC√TtargetTactual, (2.30)where the actual temperature is calculated asTactual =13Nk∑imi|vCi |2, (2.31)where k is the Boltzmann constant. Periodic temperature scaling is essen-tial in a simulation, but simply rescaling the velocities of molecules is too352.4. System Evolution“severe”.A weaker coupling method can introduce energy fluctuations in orderto simulate a canonical ensemble. A constraint thermostat, also referredto as a Gaussian thermostat [84, 86], is widely used for MD simulations.It introduces a temperature constraint into the equation of motion. Thebasic idea of this thermostat is to use a friction factor ξ to rescale velocitiesto adjust the kinetic energy of the system to the target temperature. Theequations of motion can be modified by subtracting a scaled velocity termusing the friction factor [84]aCi = Fi/mi − ξ(r,v)vCi , (2.32)and ξ is determined asξ(r,v) =∑i pi · Fi∑i |pi|2, (2.33)where pi is the momentum of particle i, and Fi is the force acting on particlei. The friction factor ξ is initially assigned to 0 and then evaluated as inEquation (2.33) at every time step. The accelerations are tuned accordingto Equation (2.32) instead of Newton’s equation of motion (Equation 2.17),and used for corrections of other properties. Note that this thermostat alsogives correct dynamical properties to linear order [86]. These calculationsare performed conveniently during the correction step of the Gear predictor-corrector algorithm.To reach the target temperature and equilibrate the system, a combi-nation of velocity scaling and Gaussian thermostat are applied during thesimulations. Usually the former method is applied every thousand steps and362.5. Comparison with Open Source Packagesthe latter one implemented in integrating the equations of motion at everytime step.2.5 Comparison with Open Source PackagesSome calculations in this work were carried out using open source packages,such as GROMACS [87] and LAMMPS [88]. As they are designed to run par-allel, their performances are especially good for large system with thousandsof molecules, and simulation runs longer than hundreds of nanoseconds. Theresults agree well with those obtained with our own programs.LAMMPS is widely used for MD simulations, and it allows us to applyvarying electric fields in the system. Our model of confined water can beset up in the LAMMPS program, and the results obtained with our owncode can be reproduced. The LJ potentials between molecular sites arecalculated straightforwardly with a cut-off at some distance, which is thesame in LAMMPS as in our program. Whereas for Coulombic interactions,the Particle-Particle Particle-Mesh method (PPPM) is used in LAMMPSinstead of the standard Ewald summation [89, 90]. In common with Ewaldsummation, the PPPM method handles the short-range and long-rangeinteractions separately. The basic idea of the PPPM algorithm is to calcu-late the short-range interactions through a direct sum, and the long-rangeinteractions from mesh points. The short-range interactions are countedthe same as the real space term in the Ewald calculation in Equation (2.9).The long-range interactions are handled by a particle-mesh technique, whichconverts charges of particles into a mesh of density values. The electrostatic372.5. Comparison with Open Source Packagespotential due to the charge distribution on the mesh is determined by solvingPoisson’s equation using Fourier transform techniques, and then the force ona charge is calculated from the mesh field by interpolation. This method isfaster than Ewald by introducing fewer mesh points than particles, and givesresults as accurate as Ewald. Water-surface potentials are also identical inall programs, where the oxygen of water molecules interact with the surface.For six-site and TIP4P/Ice water models, the LAMMPS package usesquaternions for rigid molecules [91] by computing body forces and torquesfrom atomic forces. Differing from the Gear predictor-corrector algorithm inour code, LAMMPS uses the velocity-verlet algorithm [92] for rigid molecules.LAMMPS adopts the Nose´-Hoover thermostat to probe a correct canonicalensemble efficiently. Similar to the constraint thermostat used in our pro-gram, the Nose´-Hoover method [93, 94] introduces a thermal reservoir anda friction term in the equations of motion.The results of our program are reproducible with the LAMMPS packagefor the same system. For example, the simulations in Chapter 3 wereperformed and compared for both LAMMPS and our programs. We getthe same structure of ice, the same density profile, and the same dynamicalproperties. Nucleation rates may vary for different trials of the same systemeven in the same program, so it is hard to compare quantitatively.For bulk water with uniform electric fields considered in Chapter 5,we use the programs of GROMACS. It is suitable for MD simulations oflarge system with a uniform electric field. In GROMACS, a truncated LJpotential is used for interactions between atoms and PPPM for electrostaticinteractions amongst charges. In contrast with the LAMMPS package,382.6. Data AnalysisGROMACS uses the LINCS (linear constraint solver) algorithm for rigidwater models [95]. Differing from the Gear predictor-corrector algorithm inour code, GROMACS uses the leap-frog integrator [96]. For simulations withconstant temperature and constant pressure, the GROMACS package usesthe Nose´-Hoover thermostat for constant temperature, and the Parrinello-Rahman barostat for constant pressure. [97, 98].By comparing the water and ice structures, we confirmed that GRO-MACS and our own program are consistent with each other.2.6 Data AnalysisDuring the production simulation runs, trajectories and velocities are recordedfor selected time steps. These data are read and studied through analysisprograms to obtain structural and dynamical properties. In this section,how we analyze the simulation data is discussed.2.6.1 Radial Distribution FunctionThe radial distribution function [g(r)] is a measure of the probability offinding a particle at a distance r from a reference particle, and is widelyused to describe the average structure of a system. To calculate an atom-atom g(r), a histogram Hist(r) is constructed by counting every pair ofparticles at the displacement r, and then normalized to get g(r) through[77]g(r) = Hist(r)/N4piρ/3[(r + δr)3 − r3] , (2.34)392.6. Data Analysiswhere N is the total number of atoms considered, ρ is the bulk numberdensity of the species, and δr is the radial shell thickness of 0.01 A˚. In thisthesis, the radial distribution function is usually obtained as a statisticalaverage over 1000 “independent” configurations of a production run of morethan 2 nanoseconds.As a system evolves, g(r) provides clear evidence of the structural changebetween fluid and crystal. Since liquid water and ice have different peaks,g(r) is a useful tool to track ice nucleation, but is not sensitive enough toreflect the subtle differences between cubic and hexagonal ice structures.Another algorithm for distinguish these two ice structures is introduced inthe next section.2.6.2 The CHILL Ice DetectorThe CHILL algorithm was developed by Moore et al, [75] to analyse icestructure and identify cubic and hexagonal ice. It uses the correlation oforientations among the first shell of water neighbours to classify moleculesas cubic ice, hexagonal ice, intermediate ice, or liquid water molecules. Toanalyse a specified configuration, we calculate the orientational alignmentsfor each water molecule i and its four closest neighbours ja(i, j) =l∑m=−lqlm(i)q∗lm(j)(l∑m=−lqlm(i)q∗lm(i))1/2(l∑m=−lqlm(j)q∗lm(j))1/2, (2.35)where q∗lm is the complex conjugate of qlm. qlm is defined as an instantaneouslocal order parameter, and for a particular value of l it can be calculated402.6. Data Analysisthrough m(= 2l + 1) spherical harmonics Ylm(rˆij)qlm(i) =144∑j=1Ylm(rˆij), (2.36)where rˆij is the unit vector connecting water molecule i with one of its fourclosest neighbours j. It is arbitrarily chosen that l = 3, which provides thebest resolution of ice structures and involves the least calculations [75].Water molecules can be classified according to the local structure withrespect to a(i, j) of its four neighbours. Both cubic and hexagonal ice havefour tetrahedrally coordinated neighbours in the first coordination shell andeach of the four neighbours also has a tetrahedral shell. But there is a smalldifference in symmetry properties. Molecules classified as cubic ice containonly staggered or “chair” configurations (a(i, j) < −0.8), while hexagonal icehas three staggered forms and one eclipsed, or “boat” configuration (−0.2 <a(i, j) < −0.05) [75, 99]. In this thesis, we regard both intermediate ice andliquid as liquid water instead of distinguishing them, and pay particularattention to the two ice structures.2.6.3 Orientational Order ParameterAn orientational order parameter is used to describe the orientational align-ment of the system [64]. For a given configuration, it is defined by thevectorm(z) = 1N(z)µwN(z)∑i=1µi , (2.37)412.6. Data Analysiswhere N(z) is the number of water molecules in a rectangular volumeelement centred at z, µw is the dipole moment of the particular water modelemployed, and the vector sum is over all water dipoles in the volume element.In the present calculations, the x and y dimensions of the volume elementare those of the simulation cell, and the thickness is 0.01 A˚. This orderparameter is convenient to show how dipole orientation is impacted by anelectric field. For an isotropic sample of random orientations, the meanvalue of the order parameter < m(z) >= 0; whereas for a perfectly alignedsample, < m(z) >= 1. As with g(r), the orientational order parameterprofiles shown here are averaged over 1000 configurations.In addition to the structural properties discussed above, mass densityprofiles are calculated by dividing the system into rectangular volume el-ements, and counting the number of water molecules in each rectangularelement. Plots of density profiles show the extension of ice growth. Molec-ular dynamics simulation also allows us to determine the time dependentproperties of fluids, such as the diffusion coefficient. For fluid systems,the mean square displacement < |r(t) − r(0)|2 > is a common measureof diffusion and is plotted as a function of time t [77]. The mean squaredisplacement (MSD) contains information on the molecular diffusivity. Forliquid system, MSD grows linearly with time, and the slope is the diffusioncoefficient DD = limt→∞16t < |r(t)− r(0)|2 >, (2.38)where t is the time interval, the r(t) is displacement of a molecule fromthe beginning to t, and < · · · > denote the average over all molecules422.6. Data Analysisand configurations. The mean square displacement as well as the diffusioncoefficient are useful to characterize the motional behaviour of the system.43Chapter 3Heterogeneous IceNucleation Induced bySurface Electric Fields3.1 IntroductionAs discussed in Chapter 1, heterogeneous ice nucleation is an importantphenomenon influencing many aspects of our physical environment, rang-ing from atmospheric processes to biological systems [1–3, 11, 24–26, 100].Despite its obvious importance and interest, heterogeneous ice nucleation ispoorly understood, and remains a subject of very active research. A majoradvance would be to gain an understanding of the microscopic feature(s) thatcreate a good ice nucleation catalyst. The present study is a contributiontoward this objective.One argument that is often put forward [7] invokes a good crystallo-graphic match with hexagonal ice (ice Ih) as a possible explanation ofeffective ice nuclei (IN). However, based on recent computer simulationstudies [21, 64–67], this explanation does not appear to hold when closely443.1. Introductionexamined at a microscopic level. It appears that, at least in simulationmodels, atomistically smooth, defect-free surfaces are not effective IN, evenif their crystallographic parameters do match those of ice. It is also worthnoting that a recent experimental study [68] of kaolinite particles is moreconsistent with an “active site” nucleation mechanism, than with otherpossible explanations.Another explanation that has been put forward for ice nucleation isthe possible importance of electric fields [44, 45, 47, 53, 55–58]. Computersimulations have demonstrated [47, 58, 61] that model bulk water readilyfreezes on nanosecond time scales if the sample is polarized by a strongelectric field. This is also true of confined water [62]. However, one shouldkeep in mind that in previous electrofreezing simulations all water moleculesinteract directly with the field. Thus, the field is not acting merely tocatalyze ice nucleation, but it can also influence the freezing point and theice structure obtained. Note that the ice obtained in bulk electrofreezingsimulations is ferroelectric and has the cubic ice structure (ice Ic).In previous investigations of water on kaolinite, it was observed that forsome kaolinite structures, the water near particular surfaces can be highlypolarized by strong local fields [63]. Motivated by this observation, and bythe earlier work on electrofreezing [47, 58, 61, 62], we examine field-inducedice nucleation for a simple model where water molecules experience a strongexternal field only if they are very near a surface. For a particular set ofconditions, we demonstrate that such a surface field is indeed sufficient toinduce ice nucleation.This chapter reports the results of field-induced ice nucleation, and453.2. The Model and Simulation Methodsignificantly extends the work on parameters affecting the ice nucleationprocess. We examine the important influences of temperature and electricfield parameters (strength and extent) on ice nucleation and growth. Toensure that our observations are not an artifact of a particular water modelor choice of system parameters, results are reported for two different watermodels, and for a range of system sizes and dimensions.The remainder of this chapter is divided into three parts. The modelsand simulation method are described in Section 3.2, the results are givenand discussed in Section 3.3, and our conclusions are summarized in Section3.4.3.2 The Model and Simulation MethodWe consider water confined between two infinite parallel plates that areperpendicular to the z axis, and separated by a distance Lz. The watersample occupies a simulation cell of dimensions (Lx, Ly, Lz), that is assumedto be periodically infinite in the x and y directions. The surfaces are assumedto be smooth and interact with the water through the potential given byEquation (2.1) in Chapter 2. The water-wall interaction is relatively weakand imposes little order on water near the surface [76]. This is a desirablefeature for the present work, because insofar as possible we wish to isolatethe influence of electric fields from other possible surface effects.In order to check that our results are not strongly dependent on aparticular water model, simulations were performed for two models, specif-ically, the so-called six-site model introduced by Nada et al. [72], and463.2. The Model and Simulation Methodthe TIP4P/Ice model developed by Abascal et al. [71], as described inChapter 2. Both of these water models were developed with the aim ofobtaining a better description of ice and melting. In these rigid models thewater molecules interact via site-site potentials described by combinations ofLennard-Jones (LJ) and Coulombic interactions. The TIP4P/Ice model is are-parameterization of the earlier TIP4P model [59], and the six-site modelis a recent water model [72]. Parameters for the six-site and TIP4P/Icemodels are given in Section 2.1 in Chapter 2.For the present work, it is important to know the melting points of themodels employed. For the six-site and TIP4P/Ice models the melting pointsof ice Ih at 1 atm have been estimated to be ∼ 289 K (Ref. [101]) and ∼ 270K (Refs. [71] and [102]), respectively.Mineral particles such as kaolinite that serve as IN in the physical en-vironment tend to have very rough surfaces, and one might expect localfields that vary greatly both in magnitude and direction. We consider aparticular situation where a field acts over a narrow region near a surface,and is directed parallel to that surface. This choice of geometry was initiallymotivated by our earlier observation [63] that in kaolinite trenches (partic-ular structures associated with partially cleaving basal planes of kaolinite)there tends to be a strong electric field component acting parallel to thetrench walls.In the model of field used in this chapter, the field E(z) is given thefunctional form473.2. The Model and Simulation MethodE(z) = Emax( 11 + ea(z−c))yˆ , (3.1)where Emax is the maximum magnitude of the field, z is the distance from awater oxygen atom to the surface, yˆ is a unit vector along the y axis, and aand c are parameters determining the extent and decay rate of the field asa function of z. In our model, the field is determined by the position of theoxygen atom, and is constant at that value for all point charges located inthe same water molecule. Thus the field exerts a torque but no translationalforce on a water molecule. In real physical situations, we would expect aspatially varying field to exert some force that could influence the densityand packing of water molecules near a surface, and this could vary widelydepending on the local features of a rough surface. However, since the detailsof the surface variation implied by Equation (2.15) are arbitrary, here weinclude only the polarizing aspects of the field, which we would expect tobe always present and less sensitive to details of the field gradient.All simulations are carried out under NV T conditions. In the simula-tions reported in this chapter the temperature was kept constant employinga Gaussian isokinetic thermostat [77, 86]. However, some simulations werealso carried out with a Nose´-Hoover thermostat [77], and no significantdifferences were observed. Specifically, ice nucleation and growth occurredin a similar fashion for both thermostats. A rectangular simulation cellwith Lx = Ly is used in all cases, and, as noted above, periodic boundaryconditions are applied in the x and y directions. The system is finitein z, and to account for this the water densities reported are estimated483.3. Results and Discussiontaking the cell volume to be LxLy(Lz − 2σ). The equations of motionare integrated employing a fifth-order Gear predictor-corrector algorithmwith a time step of 2 fs. The LJ interactions were cut and shifted at Lx/2to avoid discontinuities. The Coulombic interactions in the slab geometrywere treated using Ewald sums as described in Chapter 2. In the Ewaldcalculations, the parameter κ = 6.0/Lx, 1200 reciprocal lattice vectors wereincluded, and the real space terms were truncated at Lx/2. Test runsperformed with larger numbers of reciprocal lattice vectors did not showany significant differences.In all cases the sample is first equilibrated in the liquid state well abovethe melting point. The system is then cooled to the target temperature andre-equilibrated as a metastable liquid before the field is applied. We notethat even with many trials and runs of many nanoseconds, we never onceobserved ice nucleation in the absence of a field.3.3 Results and DiscussionIn this section, we demonstrate that an electric field, which acts very neara surface, can create an effective ice nucleus in models of supercooled liquidwater. We also examine the influences of the system size and cell dimensions,temperature, field strength, and field extent on ice nucleation.3.3.1 Ice Nucleation Induced by a Surface fieldIn this section, we present a particular simulation of heterogeneous icenucleation induced by a surface field, in a system that has 1200 water493.3. Results and Discussionmolecules. The density of this system is 0.96 g/cm3, and the simulationis performed at 270 K. Here we use a = 10 and c = 10 A˚ in Equation (3.1),which gives a field that falls very rapidly becoming insignificant beyond∼ 10 A˚. This field is plotted in Fig. 2.4 in Chapter 2. Larger values of cgive results similar to those obtained with c = 10 A˚, whereas smaller valuesdid not function as IN under the conditions considered here. We emphasizethat Equation (3.1) is simply a convenient functional form for a field thatdecays with distance from a surface, it is not meant to accurately representany particular physical situation.Configurational snapshots that illustrate ice nucleation and growth areshown in Fig. 3.1. At t = 0 (top panel), just before the field is applied,the system is a supercooled liquid and no order is discernible. At 2.4 ns(middle panel) after the field is applied to molecules within ∼ 10 A˚ of thesurface, ice nucleation near the surface and some growth into the bulk canbe clearly seen. At 4 ns (bottom panel), the ice has reached the vicinity ofthe opposite wall. Identical systems started from different initial conditionsshowed different time-dependent nucleation and growth behaviour, but icedid nucleate and grow in all examples (six with 1200 particles) considered.In order to understand the nature of the ice formed in our system, it isuseful to examine the dipole order parameter profile along the z axis, andthe water-water radial distribution functions.For a specific configuration, the dipole order parameter m(z) is definedin Section 2.6.3 in Chapter 2. The average components of m(z) (< mi(z) >,i = x, y, z) are plotted in Fig. 3.2. We see that in the field region, the sampleis highly polarized in the y direction, as we would expect, but beyond this503.3. Results and DiscussionFigure 3.1: Configurational snapshots of one particular simulation with 1200water molecules. The field (directed upwards along the y axis) was appliedat t = 0 (upper panel), the middle (t = 2.4 ns) and bottom (t = 4 ns) panelsshow the nucleation and growth of ice. The oxygen atoms of the watermolecules that experience the field (z . 10 A˚) are green, those outside thefield region are red, and all hydrogen atoms are black.513.3. Results and Discussionregion the polarization rapidly dies. Thus, the narrow ice layer directlynucleated by the field is ferroelectric, but the subsequent layers are not.This is also evident in the configurational snapshots shown in Fig. 3.1.Looking closely, one can see that the protons in the surface-field region arepredominately directed upwards, whereas beyond this region the protonsshow no orientational preference.−1 0 1 0  10  20  30  40<mi(z)>z(Å)Figure 3.2: The dipole order parameter profile along the z axis for the field-nucleated, frozen sample. The red, green, and blue lines are the x, y, andz components, respectively. Note that the y component shows field-inducedpolarization for z . 10 A˚.Oxygen-oxygen and hydrogen-hydrogen radial distribution functions areshown in Fig. 3.3. Results are plotted for the same sample converged underfield-nucleated, constant field (the field is uniform throughout the sample)and zero field conditions. The corresponding functions for bulk water frozen523.3. Results and Discussioninto ferroelectric cubic ice under a constant field are included for comparison.The zero field curves resemble those typically obtained for liquid water. Forthe constant field and field-nucleated systems, the oxygen-oxygen curvesare very similar, and, moreover, they closely resemble the correspondingradial distribution function of bulk cubic ice. This strongly suggests thatthe ice formed in our samples is also cubic, and we have confirmed thisby inspecting configurational snapshots along different crystal axis. Thehydrogen-hydrogen radial distribution functions for the constant field andfield-nucleated systems differ significantly, reflecting the fact that the field-nucleated ice is not ferroelectric (proton ordered) beyond the field-polarizedregion. Our observations indicate that proton disordered cubic ice, growsonto the proton ordered, ferroelectric layer, in a near seamless manner.Finally, it must be noted that in bulk water cubic ice is less thermody-namically stable than the usual hexagonal form, ice Ih. At first sight, oursimulation results would appear to be at odds with this fact. However, cubicice, that under certain conditions can anneal to form hexagonal ice, has beenexperimentally observed [103]. Also, it has been argued [104] that cubic iceis actually the preferred crystal structure for water confined to pores orto thin films (. 100 A˚ thick). Thus, our simulations suggest a plausiblescenario for heterogeneous ice nucleation. The very rough particles that actas IN in the physical environment might have some surface structures wherethe local fields are strong enough to polarize water near the surface. Indeed,given the structural variations in most real particles that serve as IN, thispossibility does not appear unlikely. Such features would serve as “activesites” for cubic ice nucleation, which could eventually convert to the more533.3. Results and Discussion 0 1 2 3 2  4  6  8g HH(r)z(Å) 0 1 2 3 4 5g OO(r)Figure 3.3: Oxygen-oxygen and hydrogen-hydrogen radial distributionfunctions. The black, red, and dark blue lines are the zero field, constantfield, and field-nucleated results, respectively. The light blue lines representferroelectric bulk cubic ice frozen under a constant field.stable hexagonal form as it grows into the bulk.3.3.2 Ice Nucleation in Systems of Varying Size andDimensionIt is important to verify that any observation of field-induced ice nucleationand growth is not unduly influenced by finite system size. Therefore, wecarried out a number of simulations for both the six-site and TIP4P/Ice543.3. Results and Discussionmodels varying the cell dimensions and the number of water moleculesinvolved. Additionally, for each system several trial runs (not less than four)beginning with different initial conditions were performed. The parametersof the systems considered, together with the number of trials, and the totallength of the runs, are summarized in Table 3.1. The number of trials thatresulted in ice nucleation for a given run length is also indicated in Table3.1. Given that we are observing ice nucleation and growth along the z axis,the surface-surface separation Lz is an important variable, and is variedbetween 40 A˚ and 80 A˚ in our simulations. It is important to test that ourqualitative observations are independent of the surface separation, becauseif Lz is too small the field could influence the freezing point itself, in whichcase we could be observing something more closely akin to electrofreezingthan to simple nucleation. The Lx = Ly dimension is also varied, and thenumber of particles in the central cell ranges from 432 to 1800. In all cases,the density is held fixed at 0.96 g/cm3.Lx = Ly, Lz (A˚) Number of Molecules Trials, Time (ns)19.6, 40 432 12(12), 32 ; 8(7), 3232.68, 40 1200 6(6), 16 ; 4(3), 1619.16, 60 648 6(6), 16 ; 4(3), 1531.9, 60 1800 4(2), 8 ; 4(4), 819.33, 80 900 9(9), 22 ; 9(9), 26Table 3.1: The systems considered in tests of possible finite size effects. Inthe third column, the first and second pairs of integers give the number oftrials and the run times for the six-site and TIP4P/Ice models, respectively.The number of trials for which ice nucleation was observed in the givensimulation time is given in brackets.In the simulations included in Table 3.1, the parameters used in Equation553.3. Results and Discussion(3.1) are a = 10 A˚−1, c = 10 A˚, and Emax = 5.0×109 V/m. As noted above,these parameters produce a field that falls rapidly to insignificance at ∼ 10A˚ from the wall, but is sufficiently strong to readily induce electrofreezing inbulk systems [47, 58], and ice nucleation in our geometry. The simulationssummarized in Table 3.1 were carried out at 270 K for the six-site model andat 260 K for TIP4P/Ice. Given the different melting points of these models,these temperatures correspond to undercoolings of ∼ 20 K and ∼ 10 Kfor the six-site and TIP4P/Ice models, respectively. The density, estimatedtaking the cell volume to be LxLy(Lz − 2σ) as discussed above, was fixedat 0.96 g/cm3, a value where bulk water has been shown [47, 58] to readilyform cubic ice for the field strength given above. This density was used forall results reported in this Chapter. An analysis of the influences of varyingthe temperature, and the field parameters is given below.For both water models, ice nucleation and growth was observed for allsystem sizes. The qualitative behaviour pattern was similar for all systemsizes and dimensions, providing evidence that finite size effects are not asignificant issue. Configurational snapshots for two examples with Lz = 60A˚ (648 particles) and 80 A˚ (900 particles) are shown in Fig. 3.4, exhibitingtwo different crystallographic planes of cubic ice. We note that in theseselected runs of 18 ns, the ice layer has grown to approximately the samethickness in both systems. However, as one would expect, the time requiredto observe ice nucleation and growth varies considerably for different trials ofthe same system. The number of trials for which ice nucleation and growthwas observed for a given run length is reported in Table 3.1. Once nucleationhas occurred the ice grows through the bulk toward the opposite surface, as563.3. Results and DiscussionFigure 3.4: Configurational snapshots after 18 ns of two selected trial runsthat show ice nucleation and growth. In one system (top), Lz = 60 A˚ (648particles) and in the other (bottom) Lz = 80 A˚ (900 particles). The oxygenatoms of water molecules that experience the field are blue, those outside thefield region are red, and all hydrogen atoms are black. The top and bottomprojections can be recognized as the (101) and (001) crystallographic planesof cubic ice, respectively.can be seen in Fig. 3.4.Ice nucleation and growth can be detected and followed by plotting thewater-water interaction energy as a function of time. Example plots foreight trial runs of the six-site model with 432 particles are shown in Fig.3.5. We observe that the interaction energy decreases a little when the fieldis turned on at 2 ns. It then oscillates about a near constant value untilice nucleation occurs, at which point the energy decreases and continues todecrease until essentially the entire sample is frozen. Note that for the eight573.3. Results and Discussiontrials shown in Fig. 3.5, the time required to observe ice nucleation variesfrom ∼ 2.5 ns to ∼ 22.5 ns.0 5 10 15 20 25time (ns)-50-48-46-44-42Energy (KJ/mol)Figure 3.5: The water-water contribution to the configurational energy forthe six-site model at 270 K. Results for eight trial runs of the 432 particlesystem are shown. The field was applied at 2 ns. Note that all eight trialsshow clear nucleation and freezing, as signalled by a rapid drop in the water-water interaction energy.The field-nucleated ice formed in our model systems is cubic, as is thecase for electrofrozen bulk ice [47, 58]. The oxygen-oxygen radial distributionfunctions, examples of which are plotted in Fig. 3.6 and compared with theresult obtained for bulk cubic ice electrofrozen under similar conditions,provide some evidence for this structure. The cubic ice structure can also583.3. Results and Discussionbe seen by inspection of configurational snapshots, as was previously donefor electrofrozen bulk samples [47, 58]. Note in particular that the (001)crystallographic plane shown in Fig. 3.4 is a clear indication of cubic ice.Additionally, we examined a number of configurations for different systemsizes employing the so-called CHILL algorithm, previously used by Mooreet al. [75] to distinguish between cubic and hexagonal ice. This confirmedthat the ice nucleated and grown in our model systems is cubic. A furtherdiscussion of cubic versus hexagonal ice is given in section 3.4, below.Although we obtain the same cubic structure, the ice obtained in thepresent field-nucleated systems differs from the bulk electrofrozen case (whereall water molecules experience the field), in that only the narrow regionthat directly experiences the field is ferroelectric. The ice grown beyondthis region is not polarized. This can be seen by calculating the dipoleorder parameter profile. In our model, the field has no x or z component,and the corresponding components of m(z) average to ∼ 0. The average ycomponent, < my(z) >, is shown in Fig. 3.7 for two trials of the 432 particlesystem, and we see that the order parameter is ∼ 1 in the field region, butfalls rapidly to ∼ 0, as the field decays at ∼ 10 A˚ from the surface. Thus,non-polarized ice nucleates and grows outward from the ferroelectric layer.For both frozen samples, density profiles along the z direction are alsoincluded in Fig. 3.7. These are obtained using the same rectangular volumeelements described above for the polarization profiles. Interestingly, thedensity profiles obtained for the different trials are not the same. Both showregular density oscillations, but the interlayer spacings are different. Thisdemonstrates that fixing the field direction and the cell geometry does not593.3. Results and Discussion 0 1 2 3 4 5 2  4  6  8g OO(r)r (Å)TIP4P/Ice 0 1 2 3 4 5g OO(r)six−siteN = 432 N = 900 N = 1200 N = 1800 N = 512 (bulk) Figure 3.6: Oxygen-oxygen radial distribution functions for the six-site (toppanel) and TIP4P/Ice (bottom panel) models. Curves for different systemsizes are shown, and compared with the bulk cubic ice result.fix the crystal axes of the ice formed, different orientations with respect tothe surface remain possible. This is further illustrated in Fig. 3.8, whereconfigurational snapshots “taken” from the same perspective are shown forboth ice structures. Clearly the crystal axes have different orientationswithin the simulation cell. Other orientations are possible, and in somecases regular oscillations are not seen in the density profile along z. In suchcases, the density profiles cannot be used to track ice nucleation and growth,603.3. Results and Discussion0 10 20 30 4000.511.522.53Density profile (g/cm3 )0 10 20 30 40z(Å)-1-0.500.51<my (z)>Figure 3.7: Density (top panel) and dipole order parameter (bottom panel)profiles for two trial runs of the six-site model with 432 particles.and one must rely on that energy and/or configurational snapshots.3.3.3 Influence of TemperatureIt is interesting to establish the temperature range over which ice nucleationand growth can be observed in simulations. Therefore, calculations werecarried out over a range of temperatures for both models. In these calcula-tions, the parameters in Equation (3.1) were again held fixed at a = 10 A˚−1,c = 10 A˚, and Emax = 5.0× 109 V/m. All simulation runs were equilibratedabove the freezing point, and cooled to the target temperature before the613.3. Results and DiscussionFigure 3.8: Configurational snapshots of the two trial runs for which densityprofiles are shown in Fig. 3.7. Both snapshots are of the configuration atthe 32nd nanosecond, and both are from the same perspective. The oxygenatoms of the water molecules that experience the field are blue, those outsidethe field region are red, and all hydrogen atoms are black. The ice crystalsclearly have different orientations in the simulation cell.field was turned on. In the six-site case, simulations were carried out withboth 432 and 1200 particles (see Table 3.1 for details) for all temperatures,and our observations were consistent for both system sizes. This is alsotrue for the TIP4P/Ice model for 260 K and higher temperatures. At lower623.3. Results and Discussiontemperatures only 432 particle systems were considered for the TIP4P/Icemodel, but given that no significant system size dependence was observed atother temperatures, or for the six-site model which was tested everywhere,we would not expect any in the TIP4P/Ice systems at low temperatures.The number of trials and the run lengths varied depending to some extenton how readily (or not) ice nucleation and growth could be seen.six-site TIP4P/IceTemp. (K) Ice Nuc. Temp. (K) Ice Nuc.230 No 230 No240 No 240 No250 Yes 245 Yes260 Yes 250 Yes270 Yes 260 Yes280 Yes 265 Yes285 No 270 Yes290 No 275 NoTable 3.2: Temperatures considered for the six-site and TIP4P/Ice models.Yes and No indicate whether ice nucleation and growth was observed or notat the particular temperature.The temperatures considered and our qualitative observations are sum-marized in Table 3.2. If “Yes” is entered in the ice nucleation (Ice Nuc.)column, then ice nucleation and growth was observed for at least one trialat that temperature. If two system sizes were considered, as was usuallythe case (see above), then this is true of both. If “No” is entered, thenice nucleation and growth was not at all observed at the temperature inquestion. Of course, it is always possible that with more and/or longersimulation runs ice nucleation might occur, but one can at least concludethat ice nucleation is much less likely outside of the temperature ranges633.3. Results and Discussionindicated in Table 3.2. 1 2 3 4 5 0  2  4  6  8g OO(r)r (Å)240 K260 K270 K280 K290 K 1 2 3 4 0  10  20  30  40MSD × 1020 (m2 )time (ps)Figure 3.9: Oxygen-oxygen radial distribution functions and mean squaredisplacement (MSD) curves for the six-site model at different temperatures.643.3. Results and DiscussionFrom Table 3.2, we see that ice nucleation and growth was observed from250 - 280 K for the six-site model, and from 245 - 270 K for TIP4P/Ice.Noting again that the normal melting points of hexagonal ice (ice Ih) havebeen estimated to be ∼ 289 K for the six-site model [101] and ∼ 270 K forTIP4P/Ice [71, 102], we would not expect to see ice nucleation above thesetemperatures, and in fact we do not. However, it should be noted that ourmodel system forms cubic ice (ice Ic) rather than the hexagonal form. Themelting points of cubic ice for the six-site and TIP4P/Ice models are notknown, but they are expected to be to be lower than those of hexagonal ice.In fact for the six-site model, cubic ice has been shown [72] to be less stablethan the hexagonal form, consistent with real water. This obviously raisesquestions as to why cubic ice prefers to grow in our model, and a furtherdiscussion of this issue is given below (Section 3.4).The fact that we did not observe ice nucleation and growth at 240 Kand lower temperatures might be puzzling at first sight, but it is in alllikelihood an artifact of the finite time (32 ns) of our simulation runs. Atlow temperatures, our model liquids become very glass-like and remain inamorphous disordered states on simulation time scales. It is possible thatmuch longer simulation runs would extend our ice observation to slightlylower temperatures, but we haven’t pursued this possibility, since it is un-likely that anything new would be learnt from such efforts. Suffice to say,that on convenient simulation times scales ice nucleation and growth canbe readily observed over the range 250-280 K for the six-site model, and245-270 K for TIP4P/Ice.Oxygen-oxygen radial distribution functions obtained for the six-site653.3. Results and Discussionmodel at different temperatures are plotted in Fig. 3.9. We note that at290 K and 240 K, the radial distribution functions do not show any ice-likestructure. At 290 K, the model is above the freezing point and remainsa stable liquid, and at 240 K we obtain a “frozen” amorphous state, asdiscussed above. Further evidence for this is provided by Fig. 3.9, wherethe molecular mean square displacements (< |r(t) − r(0)|2 >) for differenttemperatures are plotted as functions of time t. We note that diffusion inthe amorphous system at 240 K is very low, less in fact than that foundin the ice formed at higher temperatures. The radial distribution functionsand mean square displacement plots for the TIP4P/Ice model exhibit similartemperature dependences.3.3.4 Influence of Field Strength and ExtentAn electric field near a surface can create an effective ice nucleus in models ofsupercooled liquid water. To serve as an ice nucleus, the field must polarizeonly a very thin water layer (∼ 10 A˚), and the field strength required isrealistic on the relevant length scale. Our results support the idea that localelectric fields could play a major role in heterogeneous ice nucleation, par-ticularly for the very rough particles with many surface structure variations,that serve as ice nuclei in environmentally realistic situations.Clearly, an electric field is a very effective agent of ice nucleation, andit is of interest to ask how this effectiveness varies with field strength andextent. Note that by field extent we mean the distance from the surfacewhere the field maintains significant strength. In our simulations, the fieldis modelled by Equation (3.1), with Emax and c effectively determining the663.3. Results and Discussionfield strength and extent, respectively. With the parameter a fixed at 10A˚−1, the field remains nearly constant up to a distance c, and then decays,rapidly becoming insignificant. The rapid drop in field is mirrored by a rapiddrop in the dipole order parameter, as discussed above (Fig. 3.7).A number of calculations were carried out with field strengths below5 × 109 V/m, and with c values of 10 and 20 A˚. These calculations wereperformed at 270 K for the six-site model and at 260 K for TIP4P/Ice.Most simulations were done with 432 particles, but runs with systems of648 and 1200 particles (see Table 3.1 for details) were performed for somecases, and again our qualitative observations proved to be independent ofsystem size. If ice nucleation was not observed, trial runs were continued upto 32 ns.Ice Nuc.c (A˚) Emax × 10−9 (V/m) six-site TIP4P/Ice10 1.5 No No10 2.5 No Yes10 3.5 Yes Yes20 1.5 Yes Yes20 2.5 Yes YesTable 3.3: The dependence of ice nucleation on the parameters Emax andc in Equation (3.1). The simulations were at 270 K for the six-site model,and at 260 K for TIP4P/Ice. Yes and No indicate whether ice nucleationand growth was observed or not for the given field parameters.Our observations are summarized in Table 3.3. For c = 10 A˚, icenucleation was observed for Emax = 3.5 × 109 V/m and higher values forthe six-site model, and at Emax = 2.5 × 109 V/m and higher values forTIP4P/Ice. For c = 20 A˚, ice nucleation occurred for Emax = 1.5 × 109673.3. Results and DiscussionV/m and higher values for both water models. Note that we have not finetuned the nucleation thresholds by examining a finer grid in Emax, so thelimiting values might be a little lower than those quoted above. Also, thelow field limits obtained for c = 20 A˚ are very close to those we find forthe electrofreezing of bulk water, where all molecules experience the field.Thus, we would not expect the field strength threshold for ice nucleation tobe significantly reduced by allowing the field to act further than 20 A˚ fromthe surface. 1 2 3 4 5 0  2  4  6  8g OO(r)r (Å)TIP4P/ice 0 1 2 3 4 5 6six−site4.5, 10Å3.5, 10Å2.5, 10Å3.5, 20Å2.5, 20Å1.5, 20ÅFigure 3.10: Dipole order parameter profiles and oxygen-oxygen radialdistribution functions at different values of Emax and c for both the six-site and TIP4P/Ice models. The values of Emax × 10−9 V/m and of c areindicated in the figure.Radial distribution functions and corresponding polarization plots fordifferent fields are shown in Fig. 3.10. From the radial distribution functions683.3. Results and Discussionwe see that the ice structures are practically identical for the differentfields. This is not surprising, but it does serve to emphasize that thefield acts to nucleate ice, and does not appear to strongly influence thestructure obtained. The polarization plots show that in systems where icenucleation occurs, the dipole order parameter in the field region is near one.Consider for example the c = 10 A˚, Emax = 2.5 × 109 V/m case. For theseparameters the dipole order parameter is ∼ 0.9 for the six-site model whereice nucleation is not observed, and is much closer to 1 for the TIP4P/Icesystem, where nucleation does occur. The reason for the increased dipolarorder found for TIP4P/Ice is likely due to the fact that this model has alarger dipole moment (2.43 D) than the six-site model (1.89 D), and henceinteracts more strongly with the field. Given this observation, a properlypolarizable water model, which is expected [105] to have a mean liquid-state dipole moment in the range 2.6 - 3.0 D, could be expected to show icenucleation at fields lower than those we find to be effective for TIP4P/Ice.This is also likely to be true for real water.We see from Fig. 3.10, that for the same Emax, the dipole order canbe stronger for c = 20 A˚ than for c = 10 A˚. This is particularly evidentfor the six-site model with Emax = 2.5 × 109 V/m. The reason for this isthat for c = 20 A˚ more water molecules experience the field, and their col-lective interactions enhance the field’s influence increasing the polarization.Collective interactions also explain why the dipole order parameter shows adrop very near the surface, where the water molecules have fewer adjacentneighbours with which to interact.693.4. Summary and Conclusions3.4 Summary and ConclusionsWe have examined field-induced ice nucleation for a model system whereonly water molecules very near a surface experience an external field. Twowater models, several sets of field parameters, and a range of temperaturesare considered. Additionally, we carefully check that system size and/or celldimensions do not have a significant influence on the simulation results. Ourmain findings are as follows.Field-induced ice nucleation and growth readily occurs for both the six-site and TIP4P/Ice water models. Although there are some quantitativedifferences due to differing model parameters and physical properties, thequalitative behaviour of both models is very similar. Also, our qualitativeobservations do not exhibit any significant dependence on the system size,or on the dimensions of the cell employed. The time required for ice tonucleate after application of the field varies from a few nanoseconds to a fewtens of nanoseconds for trial runs that differ only in their initial conditions.We demonstrate that although cubic ice is nucleated and grows in all cases,fixing the cell geometry and the field direction does not fix the crystal axesof the ice formed. The axes of the ice crystals formed in different trial runscan take on different orientations in identical simulation cells.Simulations were carried out for a range of temperatures with fixedfield parameters, and ice nucleation and growth was observed over thetemperature range 250-280 K for the six-site model, and over 245-270 K forTIP4P/Ice. Recalling that the normal freezing points are estimated to be∼ 289 for the six-site model, and ∼ 270 for TIP4P/Ice, these temperature703.4. Summary and Conclusionsranges are physically reasonable. The TIP4P/Ice model does appear toexperience ice nucleation closer to its freezing point than the six-site case,but given the uncertainties in the freezing point estimates, and our fivedegree temperature grid, too much significance should not be attached tothis apparent difference. On the low temperature side, both models becomevery glassy and ice nucleation is not observed on simulation time scales.In our model, the electric field is governed by two, physically importantparameters. These are the field strength Emax and the parameter c, whichis essentially the distance from the surface where the field strength remainssignificant. We find that if c = 10 A˚, then the six-site and TIP4P/Icemodels show ice nucleation for Emax & 3.5 × 109 V/m, and Emax & 2.5 ×109 V/m, respectively. If c = 20 A˚, both models shows ice nucleation ifEmax & 1.5 × 109 V/m. Ice nucleation is only observed if the dipole orderparameter in the surface region is very nearly one. The dipole moment ofTIP4P/Ice (2.43 D) is significantly larger than that of the six-site model(1.89 D), and hence the TIP4P/Ice molecules interact more strongly withthe field. This likely explains the stronger polarization that allows ice tonucleate at a lower field strength for TIP4P/Ice when c = 10 A˚. Withc = 20 A˚, more molecules experience the field and their collective interactionleads to stronger polarization, and ice nucleates at a lower field strengthfor both models. We note that increasing c further will not enhance icenucleation. This is because the minimum field strength for which we observebulk electrofreezing (all molecules experience the field) is essentially thevalue we obtain for c = 20 A˚. Also, it is worth remarking that a properlypolarizable water model (and perhaps real water) would likely nucleate ice at713.4. Summary and Conclusionsstill lower field strengths for c = 10 A˚. The mean liquid-state dipole momentof polarizable water models (2.6-3.0 D) is larger than that of TIP4P/Ice,and hence we would expect stronger polarization near the surface, and icenucleation at a lower field strength.Finally, we note that in our model system the ice formed is cubic, just asit is in electrofrozen bulk systems. The only difference is that the bulksamples are globally ferroelectric, whereas in the present case only thesurface layer that directly experiences the field is ferroelectric, the ice thatgrows beyond this region is dipole disordered. This raises questions as towhy cubic ice rather than the more stable hexagonal form grows in oursimulations. One obvious factor is that the field-induced, ferroelectric icelayer formed near the surface has a cubic structure, and we might expect thisto give a kinetic bias to the structure of the dipole disordered ice that growsbeyond the surface layer. However, recently reported experimental and sim-ulations studies suggest that the explanation might be more complex. Cubicice that slowly converts to hexagonal ice has been reported in nucleationexperiments [103], and it has since been convincingly argued [106] that icenucleated in supercooled water actually consists of both cubic and hexagonalforms in randomly stacked layers. Similar mixed cubic/hexagonal structureshave been observed in MD simulations employing the mW model [107–109].However, only cubic ice was found [110] in metadynamics simulations usingTIP4P water, so the simulation results may have some dependence on themodel and/or the conditions employed. We note that it has long beensuggested that cubic ice might have a special role in ice nucleation [111],and it has been argued that cubic ice may be more stable than hexagonal723.4. Summary and Conclusionsice in finite systems such as thin films [104]. In any case, both watermodels we consider form cubic ice in the field nucleation geometry employed.We speculate that as the cubic ice layers extend further and further fromthe surface into the bulk, the structure would eventually convert to thehexagonal form. However, it is possible that both the sample size and thetime required to observe such behaviour lie well outside what is presentlypossible with computer simulations. The issue of cubic versus hexagonal iceis further discussed in the following chapter.73Chapter 4Ice Nucleation by ElectricSurface Fields of VaryingRange and Geometry4.1 IntroductionAs discussed in previous chapters, heterogeneous ice nucleation is an im-portant process in many natural settings [3, 11, 24, 100]. An importantexample from atmospheric science is ice clouds formation, facilitated byaerosol particles, which significantly influence climate [4, 5, 7]. There aremany other examples of atmospheric [112] and biological [24, 100] systems,where heterogeneous ice nucleation is an important factor. Despite the im-portance of heterogeneous ice nucleation, it is fair to say that the microscopicmechanism(s) of nucleation, and the particular features required of good icenuclei remain poorly understood. It is often argued [7] that a good crystal-lographic match with the ice lattice, as is the case for some ice nucleatingmaterials, is an important feature, but recent computer simulation studiesof model systems have been unable to provide any evidence in favour of744.1. Introductionthis explanation [21, 64–67]. Moreover, since very varied substances rangingfrom mineral particles to proteins and bacteria can nucleate ice, a generalexplanation based on matching crystal structures seems unlikely.It is difficult to obtain direct experimental information about early stageice nucleation, due at least in part to the small length scales and shorttimes involved. An additional complication is that the particles which areimportant ice nuclei in realistic situations, such as the mineral particlesrelevant to atmospheric science, nearly always have very rough surfaces witha variety of structural features that could possibly serve as “active sites”for ice nucleation [68]. From this perspective, controlled experiments withparticles that are microscopically well characterized would likely providevaluable insight, but at present there are few if any such studies. Giventhis situation, computer simulations of model systems offer another routeto information and understanding at the molecular level, which is rapidlybecoming more viable as computational power increases.On the simulation side, there have been a number of interesting studiesof homogeneous ice nucleation employing both direct molecular dynamics(MD) methods [107, 113], and other indirect “biased sampling” techniques[99, 110, 114]. For some questions biased sampling methods are preferred,providing an efficient means of sampling rare events that may not be ac-cessible by any other simulation approach. Biased sampling methods areexcellent for obtaining thermodynamic properties such as free energy differ-ences that are independent of path, but because they require the a prioriselection of a “reaction coordinate”, such methods may reveal little aboutthe actual nucleation process. Direct MD simulation has the advantage754.1. Introductionthat ice nucleation and growth can follow the physically relevant pathway,but the disadvantage that for realistic atomistic water models homogeneousice nucleation events are too rare for reliable study on currently feasiblesimulation timescales. Direct MD simulations have proven extremely usefulfor the more coarse grained monatomic water (mW) model, developed andextensively studied by Molinero and coworkers [75]. A great deal of insightinto the nature of homogeneous ice nucleation has been obtained from thiswork [107, 113].A different situation might be expected for heterogeneous ice nucleation,where a good ice nucleation catalyst could dramatically speed up the processsuch that direct MD studies of ice nucleation and growth become possiblefor realistic water models. Perhaps the most favourable example currentlyknown is ice nucleation by external electric fields that can apply over theentire sample [47, 61, 62], or merely over a thin region near a surface[115, 116]. The case where the field applies throughout the sample is moreproperly termed electrofreezing because in that case the field influences thefreezing point itself, whereas a local surface field acts only as a catalystthat facilitates the nucleation process. In such simulations, ice nucleationinduced by an applied field is a robust phenomenon not strongly dependenton simulation variables, such as sample size. We note that while therehave been experimental reports of field influenced ice nucleation [55–57], adefinitive demonstration with well defined and measured fields is still lacking.One inhibiting problem is that the fields required (at least on simulationtimescales) are orders of magnitude larger than those that cause dielectricbreakdown in macroscopic water samples, and it is difficult to create con-764.1. Introductiontrolled fields on the small length scales required to avoid breakdown [44].Direct MD simulations employing atomistic water models have also beenused by Jungwirth and coworkers to investigate ice nucleation at free watersurfaces [70, 117], as well as at water-pentanol and water-pentanoic acidinterfaces [118]. For rather narrow simulation cells (∼ 1.35 nm) repeatedperiodically in directions perpendicular to the interface, ice nucleation isobserved in these systems, but apparently freezing is not observed for widerperiodically repeated simulation cells [118]. Recently, somewhat similarstudies of water-kaolinite interfaces were reported by Cox et al. [21], wherethey observe an interesting influence of kaolinite on ice nucleation. However,as clearly noted by the authors, the observations again strongly depend onsample size, so the status and significance of the results remains somewhatunclear.In this chapter we further explore ice nucleation by surface fields. Inearlier work [115, 116] discussed in Chapter 3, we considered ice nucleationby fields that were finite perpendicular to the surface, but acted over theentire surface area. Here we examine fields that act only over surface“bands” of different size and geometry. We show that field bands of variedshape can function as very effective ice nuclei, provided that a certain sizethreshold is exceeded. By varying the field bands, we also learn more aboutthe mechanism of ice nucleation and growth. One important observationconsistent for all systems considered is that essentially all ice growth occursat (111) planes of cubic ice. By varying the field geometry and proportionof the surface it covers, we gain insight into how field nucleation influencesthe mix of hexagonal and cubic ice observed under different simulation774.2. Models and Methodsconditions.The remainder of this chapter is divided into three parts. The modelsand simulation method are described in Section 4.2, results are given anddiscussed in Section 4.3, and our conclusions are summarized in Section Models and MethodsThe so-called six-site [72] and TIP4P/Ice [71] water models are both con-venient for simulation studies of ice and freezing. In earlier investigations[115, 116] of field-induced ice nucleation discussed in Chapter 3, we haveshown that the physical behaviour of both models is qualitatively similar.Therefore, in the present chapter we report results only for the six-sitemodel, but we have verified that no significant differences occur in the be-haviour of TIP4P/Ice. The details of the six-site water model is described inSection 2.1, and the normal freezing point of the six-site model is estimated[101] to be ∼ 289 K.As in our previous work [115, 116], the water molecules are confinedbetween two parallel surfaces as sketched in Fig. 2.2 in Chapter 2. Thesurfaces are separated by the distance Lz, which together with Lx and Lydefines the dimensions of the rectangular simulation cell. The system isperiodically infinite in the x and y directions. The surface-water interactionsare given in Equation (2.1) by Chapter 2. These parameters correspond toa “neutral” surface, which imposes little order in the neighbouring watermolecules.In our slab geometry (Fig. 2.3 in Chapter 2), water molecules near the784.2. Models and Methodsleft-hand surface experience an electric field directed along the x-axis thatis of finite range in z, and possibly in y. The spatially varying field is givenbyE(y, z) = Emax( 11 + f(y, z))xˆ , (4.1)where f(y, z) is very small in the region of the field, but grows exponentiallyfast beyond the designated boundary. Note that in this chapter the fieldis taken to be in the x direction, and in Chapter 2 (Fig. 2.3) the fieldis illustrated in y direction. Equation (4.1) gives a field which is nearlyconstant at Emax within a defined volume, but decays rapidly to zero outsidethis volume. We consider field bands (finite range in z and y) with crosssectional areas that are rectangular, triangular, and semicircular in shape,as shown in Fig. 2.3. Appropriate choices of f(y, z) control the shape andextent of the field bands.In previous simulations [115, 116] (Chapter 3), we employed our ownMD program, but in the present calculations we use the LAMMPS MDpackage [88], which is a more efficient parallel code, allowing us to simulatelarger systems for longer times. Periodic boundary conditions (PBC) areformally employed in all three directions, but to accommodate the slabgeometry empty regions are placed outside the confining surfaces locatedon the z axis [76]. The Particle-Particle Particle-Mesh (PPPM) method [90]is used to take account of the long-range Coulombic interactions. The short-range, water-water and water-wall interactions are truncated at 0.8 nm. Allsimulations are performed in the NV T ensemble. The equations of motion794.2. Models and Methodsare integrated applying the velocity Verlet algorithm [77] (with a 0.75 fstimestep) for the rigid six-site water model. The Nose´-Hoover thermostat[93, 94] with a relaxation time of 0.037 ps is used to control the temperature.In all simulations, liquid water is first equilibrated at 350 K, which is wellabove the freezing point (∼ 289 K) of the six-site water model. The samplesare then cooled to 270 K, and relaxed at that temperature for 300 ps beforethe field is switched on. This ensures that all samples are initially in thesame relaxed liquid state before being exposed to the field. The simulationswere then continued until ice nucleation and growth was sufficiently wellestablished for our analysis. The run times (after application of the field)required for this ranged from 62 to 130 ns, depending on the system, andthe particular simulation run.In the present simulations, the applied field has Emax = 3.0× 109 V/m.This is smaller than the value (5.0 × 109 V/m) used in previous work [47,115, 116], but is sufficiently large to nucleate ice in all system geometriesconsidered here. We note that in our earlier simulations [115, 116] (Chapter3), the field experienced by all charges within the same molecule was keptstrictly constant. This ensured that the field could impose a rotationaltorque, but no translational force on a molecule. It is difficult to apply thisconstraint in the LAMMPS code, so in the present simulations the field canvary a little across a molecule. However, by comparing with strictly constantfield results for some test cases, we determine that this has no significantinfluence on ice nucleation.We use the CHILL algorithm developed by Moore et al. [75] to analyzethe ice structure. This algorithm allows one to classify water molecules804.3. Results and Discussionas belonging to cubic, hexagonal, and “intermediate” ice structures, or toliquid water. This algorithm is discussed in details in Chapter 2. Here we areinterested in identifying regions of cubic and hexagonal ice, and therefore inthe analysis given below we classify water molecules as cubic, hexagonal, orliquid, with both liquid and intermediate molecules included in the “liquid”category.4.3 Results and Discussion4.3.1 Full Surface FieldBefore discussing results for fields acting over different portions of the surfacearea, it is useful to briefly revisit the full surface field case considered earlier[115, 116] (Chapter 3). In that model, a parallel field acts over the entiresurface area (periodic in x and y), but decays rapidly with distance fromthe surface z. This situation corresponds to Equation (4.1) with f(y, z) =f(z) = ea(z−c), where the parameters a and c control the decay of the fieldwith distance from the surface. Previously, we showed that a thin surfacefield that decays essentially to zero at 0.1 nm from the surface nucleates icevery effectively. Ferroelectric cubic ice rapidly forms a surface layer, andunpolarized ice grows outward from this nucleus. However, in our earlierwork we did not determine the mechanism of ice growth.Here, we address this question by carefully observing the pattern of icegrowth in a number of simulations (20 in all) employing larger systems (≥2400 water molecules) with different initial conditions, and in some examplesdifferent cell dimensions. The same growth mechanism is observed in all814.3. Results and Discussionsimulations, and is illustrated by the configurational snapshots given in Fig.4.1. For the particular case shown, Lx = Ly = 3.06 nm, Lz = 8.5 nm, andthe field parameters are a = 0.1 nm−1, c = 0.1 nm, and Emax = 3.0 × 109Vm−1. Note that in order to give a clearer picture of the behaviour of theperiodic system, Fig. 4.1 includes the central simulation cell (outlined inblack), together with the above and below periodic images. Further note,that the perspective shown is directly into the (110) plane of cubic ice.Our observations illustrated in Fig. 4.1 are as follows. Initially (∼ 10 ns),an ice sheet ∼ 0.1 nm thick forms in the field region. However, beyond thissheet the ice does not grow layer by layer parallel to the wall. Rather, thegrowing ice layers can take different orientations, which in general are notparallel to the surface. Based on careful examination of different simulations,we determine that ice growth occurs at the (111) plane of cubic ice. This isthe case in all of our simulations, independent of the particular orientationtaken by the ice crystal with respect to the simulation cell. This is especiallyeasy to see in the snapshots shown in Fig. 4.1, where the dashed greenlines (stage 1) mark the ice-water interface at the (111) crystal plane. Icegrows layer by layer at this interface. We remark that Takahashi [111]has previously pointed out that the (111) plane of cubic ice gives the mostfavourable ice-water interfacial energy, which is consistent with our empiricalobservation. Also, as described below, the (111) plane growth mechanism isobserved in all of our simulations including those where nucleation is inducedwith fields of different geometry covering different surface areas.Another aspect of our earlier simulations that requires explanation is thefact that we observed only cubic ice in our samples [115, 116], whereas based824.3. Results and Discussion1224Figure 4.1: Configurational snapshots showing ice nucleation and growth fora 2400 molecule system with a full surface field. The snapshots show the(110) plane of cubic ice at 20 ns (stage 1), 50 ns (stage 2), 60 ns (stage 3), and70 ns (stage 4). The water oxygen atoms in the field region are dark blue,those associated with cubic ice or liquid are red, and those in hexagonal icelayers are light blue. All hydrogen atoms are black. The rectangles outlinethe central cell. The dashed green lines shown in stage 1 indicate growth at(111) planes of cubic ice.834.3. Results and Discussionon other simulation [73, 75, 99, 119] and experimental [103, 106] studies ofice growth under similar conditions, some mixture of cubic and hexagonalice layers is expected. As discussed by others, cubic ice (ice Ic) is onlymarginally less stable than hexagonal ice (ice Ih), and, moreover, the basalplane of hexagonal ice and the (111) plane of cubic ice are identical [73, 74].Therefore, cubic and hexagonal ice can be “joined” at this surface withonly a small unfavourable interfacial energy. Despite this, and in contrastwith other simulations of ice growth [73, 99], in over 40 simulations of icenucleation by full surface fields, we observed a hexagonal ice layer only once.In all other cases only cubic ice was formed. However, in simulations wherethe field acts only over a portion of the surface area (see discussion below)we routinely observe both cubic and hexagonal ice layers in the growingcrystal.At first sight this discrepancy appears rather surprising, but we nowbelieve that it is an artifact of finite system size and PBC. The singleexample where we do observe a layer of hexagonal ice under full surfacefield conditions is shown in Fig. 4.1 (see stages 3 and 4). We note thatunder PBC all ice layers must terminate in the field region near the surface.This is an unfavourable situation for hexagonal ice because near the surfacethe field strongly favours cubic ice. We note that in the example shown inFig. 4.1 the hexagonal ice layer converts to cubic near the surface regiongiving rise to a structural defect. It seems highly likely that the mismatchof hexagonal ice with the cubic surface layer strongly inhibits the formationof hexagonal layers in our finite periodic samples. In systems with partialsurface fields, the field does not extend over the entire surface area, and844.3. Results and Discussionhexagonal layers can continue unaltered contacting the surface outside thefield region, which is in fact what we observe (see below).4.3.2 Partial Surface FieldsAs noted above, the main purpose of the present chapter is to determineif, and to what extent, ice nucleation is influenced by the shape and size ofthe surface field. To that end, we have performed a number of simulationswhere the field acts in rectangular, triangular, and semicircular bands asillustrated in Fig. 2.3. Details of the simulations carried out are summarizedin Table 4.1. In simulations 1-10 (Table 4.1), 3900 molecules were used withLx = 2.67 nm, Ly = 10 nm, and Lz = 5 nm, giving a density of ∼ 0.96g/cm3. Simulation 11 involves a larger system (6000 molecules) and wascarried out to investigate any influence of the surface-surface separation.For this case, the Lx and Ly dimensions (2.68 and 10 nm, respectively) areat or near those of the smaller samples, but Lz = 7 nm, again giving adensity ∼ 0.96 g/cm3. Additional tests were carried out for systems withLy increased to 15 nm (5760 water molecules) to ensure that we are indeedseeing the effect of a single field band, uninfluenced by periodic repetition inthe y dimension. The details of these results are not given because they areessentially the same as those obtained with Ly = 10 nm, providing evidencethat our observations are not significantly influenced by the periodicity iny.The dimensions (y, z) of the three field geometries considered are givenin the table. For rectangles, (y, z) represents (length, width), for triangles(base, height), and for semicircles the radii are given. The functions f(y, z)854.3. Results and Discussionsimulation time (ns) cubic hexagonal liquid1, rectangular (0.65, 0.15) 62 2354 609 9372, rectangular (0.55, 0.15) 90 2236 588 10763, rectangular (0.45, 0.15) 140 1222 28 26504, rectangular (0.45, 0.15) 90 2941 140 8195, rectangular (0.35, 0.15) 95 2087 820 9936, rectangular (0.35, 0.15) 130 2486 488 9267, triangular (0.75, 0.19) 70 2093 841 9668, triangular (0.70, 0.18) 75 1905 881 11149, semicircular (0.25) 90 2081 943 87610, semicircular (0.23) 95 2139 783 97811, rectangular (0.45, 0.15) 85 3656 1072 1272Table 4.1: Summary of simulations carried out with partial surface fields.For simulations 1-10, Lx, Ly, and Lz are 2.67, 10, and 5 nm, respectively,for simulation 11, 2.68, 10, and 7 nm, respectively. In all simulations thedensity is ∼ 0.96 g/cm3. The numbers given in brackets in column one arethe (y, z) dimensions (in nm) of the field bands (Fig. 2.3). The times givenin this table are the lengths of the simulations. For rectangles these are thelength and width, for triangles the base and height, and for semicircles theradius. The numbers of cubic, hexagonal, and liquid water molecules in theconfiguration at the final timestep are given in the last three columns.used in Equation (4.1) are selected such that the field strengths becomenegligible (i.e., much lower than that needed to nucleate ice) outside of thefield dimensions given in Table 4.1.From the results summarized in Table 4.1, we see immediately that fieldbands of all three shapes can serve as effective ice nuclei. Most trials wereperformed with the rectangular geometry, and within our simulation timesice nucleation and growth was observed for (y, z) ≡ (0.35 nm, 0.15 nm) (5out of 12 trials showed ice nucleation within 50 ns), but not for smallersizes. Simulations for (y, z) ≡ (0.29 nm, 0.15 nm) show fluctuating ice-liquidstructure in the field region, but ice growth beyond the field region was864.3. Results and Discussionnot observed. It appears that surface field bands . 0.35 nm wide are notsufficient to generate ice nuclei that live long enough to grow.Configurational snapshots showing the pattern of ice nucleation andgrowth for rectangular field bands are shown in Fig. 4.2 (simulation 2 inTable 4.1) and Fig. 4.3 (simulation 11 in Table 4.1), and for a semicircularfield band in Fig. 4.4 (simulation 10 in Table 4.1). For triangular field bandsconfigurational snapshots (not shown) indicate a pattern of ice nucleationand growth much like that of the other geometries. By and large, the icenucleation and growth process for field bands is similar to that observed forfull surface fields, and does not show any strong dependence on the shapeand size of the field band, provided that it is over a size “threshold”, asnoted above. Specifically, a block of polarized cubic ice nucleates in theregion of the field band, and unpolarized ice then grows into the bulk. Icegrowth at the (111) plane can be easily seen for the rectangular exampleshown in Fig. 4.2 (note the growth regions marked with dashed green linesin the 15 ns snapshot).Despite the similarities, systems with field bands differ in one significantway from the full surface field case. As noted above, with a full surface field itis very unusual to observe any hexagonal ice layers mixed in with the cubicstructure. This is surprising because cubic and hexagonal ice layers canbe joined with very little interfacial energy penalty, and more importantlymixed ice structures are observed in other simulation and experimentalstudies of ice nucleation and growth. As discussed above, we suspect thatthe near absence of hexagonal layers in our full surface field simulations isdue to finite system size and PBC. Under PBC ice layers terminate in the874.3. Results and DiscussionFigure 4.2: Configurational snapshots showing ice nucleation and growth fora rectangular field band (simulation 2 of Table 4.1). From left to right thesnapshots correspond to 15, and 90 ns. The atoms are coloured as in Fig.4.1. The dashed green lines in the snapshot at 15 ns indicate ice growth at(111) faces of cubic ice. Note that the hexagonal ice layers meet the surfaceoutside the field region.884.3. Results and DiscussionZYXFigure 4.3: Configurational snapshot (at 85 ns) showing ice nucleated by arectangular field band in a larger system (simulation 11 of Table 4.1). Theatoms are coloured as in Fig. 4.1. Note that in the larger system hexagonalice layers can form parallel to the surface.894.3. Results and Discussionpolarized surface ice layer, which strongly disfavours the hexagonal form.The simulations with field bands provide strong evidence that this ex-planation is correct. In nearly all simulations with field bands (Table 4.1)we find significant numbers of hexagonal ice molecules, and hexagonal icelayers can be easily identified in the configurational snapshots (Figs. 4.2,4.3, and 4.4). Moreover, we note that the hexagonal ice layers are locatedin the simulation cell such as to avoid contact with the band of polarizedcubic ice. For the smaller systems (Figs. 4.2 and 4.4), the hexagonal layerstend to be oriented across the cell from surface to surface, contacting theleft-hand surface outside the field region. For the larger system (Fig. 4.3),the hexagonal ice layers can take other orientations in the simulation cell,as we would expect in a truly infinite sample. These observations confirmour suspicion that it is the polarized surface ice layer which acts againstthe formation of hexagonal ice in our full surface field simulations. Ofcourse with sufficiently large samples, the surface mismatch and associateddefects will become thermodynamically unimportant, and we would expecthexagonal ice structures even with full surface field nucleation.We note also that in our simulations the ratio of hexagonal to cubic iceis usually less than the 1:2 ratio reported by Moore et al. [75] (see numbersgiven in Table 4.1). This difference is explained by the fact that in oursimulations all ice formed in the field band is cubic, and, furthermore, theblock of cubic ice at the surface can inhibit hexagonal ice formation in otherparts of the simulation cell.One other point worth mentioning is that the amount of hexagonal iceformed can vary considerably for simulations that differ only in their initial904.3. Results and DiscussionFigure 4.4: Configurational snapshots (at 90 ns) of ice nucleated with asemicircular field band (simulation 10 of Table 4.1). The (001) (left panel)and (110) (right panel) planes of cubic ice are shown. The atoms are colouredas in Fig. 4.1. Note that the hexagonal ice layers meet the surface outsidethe field region.914.4. Summary and Conclusionsconditions. We have noticed that in simulations where the ice grows quickly,rapidly filling the cell, there tends to be less hexagonal ice than for caseswhere the ice grows more slowly. This might indicate that cubic ice hassome kinetic advantage over hexagonal in the growth process. We also notethat once the simulation cell has completely frozen, the ratio of hexagonalto cubic ice shows no change on the timescale of our simulations.4.4 Summary and ConclusionsIn this chapter we use MD simulations to examine ice nucleation by partialsurface fields. Results are reported for the six-site water model [72] under-cooled by ∼ 20 K. For comparison some simulations were also carried outwith the TIP4P/ice water model [71], and the results obtained were notsignificantly different. We show that field bands exceeding a minimum size(& 0.15 nm thick, and & 0.35 nm wide for the rectangular case) catalyze icenucleation just as efficiently as full surface fields. Additionally, the shape ofthe field band does not appear to be particularly significant, with rectan-gular, triangular, and semicircular cross sectional areas giving qualitativelysimilar results. These observations further support the suggestion that localelectric fields might be a significant factor in some ice nuclei.Simulations with fields of different geometry and dimension also teachus more about the ice nucleation and growth process. For all systemsconsidered, based on careful observation of the pattern of crystal growth, wedetermine that nearly all growth occurs at (111) planes of cubic ice. Thisis despite the fact that the crystal axes of the growing ice can take different924.4. Summary and Conclusionsorientations in the simulation cell, varying with different initial conditions.Our simulations clearly show that growth at the (111) plane of cubic ice isstrongly favoured, consistent with a suggestion of Takahashi [111] based onconsideration of interfacial energies.It is also worth remarking that with partial surface fields we nearly al-ways observe hexagonal ice layers mixed into the cubic structure, consistentwith other simulations [73, 75, 99, 119], and with experiment [103, 106]. Thiscontrasts with the full surface field case, where hexagonal ice is almost neverobserved in our simulations [115, 116]. We show that in all likelihood thisdifference is a simulation artifact related to finite system size and periodicboundary conditions.93Chapter 5UnderstandingElectrofreezing in WaterSimulations5.1 IntroductionPure water freezes when an ice cluster grows faster than it melts. Sanz et al.[120] estimate that it would take more than 13 billion years to nucleate ice in1018 tonnes of water, at -20◦C. Therefore, a nucleating agent must usually bepresent for water to freeze. Diverse materials nucleate ice [1, 3, 10, 20, 121],but the properties responsible for ice nucleation remain unknown, or at bestpoorly understood.As we demonstrated in earlier chapters, partial electric fields can induceice nucleation efficiently. Model water freezes rapidly (within nanoseconds)when exposed to a uniform electric field on the order of 1 V/nm (i.e., ×V/m are used elsewhere), as demonstrated by molecular dynamics simula-tions [47, 58, 62]. Experiments also showed that local electric fields are apossible significant feature of many ice-nucleating materials [55–57]. Real945.1. Introductionwater undergoes dielectric breakdown if exposed to such strong, uniformfields, making it impossible to confirm simulation results experimentally[46]. However, real water can withstand large local fields produced by ions orcharge groups, and simulations indicate that local fields can also nucleate ice[115]. Some experimental studies suggest that weaker electric fields promotefreezing [55, 57, 122, 123], while others find that particular fields can inhibitfreezing [57], or have no effect [44]. Simulations can only probe, at most,microsecond timescales; if weak fields promote freezing, they may requiresignificantly more time to generate ice nuclei.The factors that enhance freezing in polarized water are unclear. Ingeneral, homogeneous freezing occurs when a critical ice nucleus forms inpure bulk water. According to classical nucleation theory, the size of thecritical nucleus decreases as either the temperature and/or the solid-liquidsurface tension decreases. Assuming classical nucleation theory accuratelydescribes ice nucleation, electric fields may promote freezing by increasingthe melting temperature, and/or decreasing the surface tension between iceand water.In this work, we find that the melting temperature of the six-site watermodel [101] increases when interacting with a uniform electric field. Inpresence of a field, polarized cubic ice coexists with polarized liquid waterat temperatures well above the normal melting point of the six-site modelat zero field. In cubic ice, the water dipoles are essentially perfectly alignedby the electric field, whereas only partial polarization is observed in thecoexisting liquid. Thus, favourable interactions with the field help stabilizeice at a higher temperature. Our simulations show that the spatial structure955.2. Simulation and Analysis Methodand tetrahedral order of liquid water change little with partial polarization;however, structural changes do occur as liquid water dipoles approach nearperfect alignment with the field. Even when strongly polarized, watercan be supercooled to ∼ 40 K below the field-dependent melting point,which is comparable with the degree of supercooling commonly observedexperimentally at zero field. This indicates that the barrier to homogeneousice nucleation remains high in the presence of an electric field. Our resultssuggest that local fields generated by ice-nucleating materials may increasethe effective melting temperature of the locally polarized water, therebyreducing the size of the critical nucleus.The remainder of this chapter is divided into three parts. The methodsof simulation and analysis are described in Section 5.2, results are discussedin Section 5.3, and our conclusions are summarized in Section Simulation and Analysis MethodMolecular dynamics simulations of the six-site water model [72] were per-formed using the GROMACS 4.5.4 simulation package [87]. As discussedbefore, the six-site water model has a high melting temperature [101] (∼289 K) and is designed to study water-ice systems near the melting point.In general, structural and thermodynamic properties of ice and water nearthe melting point are well reproduced by this water model [72].As a rigid model, six-site water molecules interact via site-site potentialsthrough both Lennard-Jones (LJ) and Coulombic interactions, as describedin Chapter 2. Short-range interactions were truncated at 0.8 or 1.0 nm;965.2. Simulation and Analysis Methodresults were qualitatively insensitive to the truncation length used. ThePPPM method was used to calculate long-range electrostatic interactions.Simulations were performed with N = 1000, 8000, or 32000 water molecules,at constant temperature (T ) and under constant pressure (P ). Constanttemperature was maintained using the Nose´-Hoover thermostat [93, 94] andthe Parrinello-Rahman barostat was used to maintain a constant pressure of1 bar [97, 98]. The equations of motion were integrated using the leap-frogalgorithm with a 2 fs time step [96].Changes in free energy and entropy due to a uniform electric field werecalculated using thermodynamic integration. In NPT simulations, the Gibbsfree energy difference (∆G) was obtained through∆G = Gfield −Gno field =∫ 10〈∂H(λ)∂λ〉λdλ, (5.1)where H is the enthalpy, the angular brackets indicate an ensemble average,and λ controls the strength of the interaction between water molecules andthe field. Ufield is calculated by summing over all charge sites in the systemUfield(λ) = −λ5N∑i=1qi(ri · E), (5.2)where E is the electric field, qi represents the i-th charge, ri is the positionvector of charge i, and 5N is the total number of charges in the system.In the presence of an electric field, H=Uconf+Ufield+PV , where Uconf inthe configurational energy due to water-water interactions, and Ufield is theenergy due to molecular charge sites interacting with the external field.975.2. Simulation and Analysis MethodSince Uconf does not depend on λ, and ∂(PV )/∂λ is two to four orders ofmagnitude smaller than ∂Ufield/∂λ, Equation (5.1), can be rephrased as∆G ≈∫ 10〈∂Ufield(λ)∂λ〉λdλ. (5.3)The entropy change due to interaction with an electric field can then beestimated,Sfield − Sno field = ∆S = (∆H −∆G)/T, (5.4)where ∆H = Hfield−Hno field. This method is valid until water freezes, andthe thermodynamic functions decrease discontinuously with ice formation.The average dipole order parameter, or average polarization, along thedirection of the applied field (x) can be calculated through< mx >=1NconfNconf∑i1NµwN∑i=1µx,i , (5.5)where Nconf ≥ 100 is the number of independent configurations included inthe average, µw is the dipole moment of the six-site water model, and µx,iis the x-component of the dipole vector of the i-th water molecule.Polarization might be expected to profoundly affect water’s structuralproperties, and as a consequence, its thermodynamic properties. To under-stand structural changes that occur in water due to polarization, we compareradial distribution functions, g(r), structure factors, S(k), and the localtetrahedral order parameter [124], q, in polarized and unpolarized water.985.2. Simulation and Analysis MethodThe local tetrahedral order parameter [124]q = 1− 383∑i=14∑j=i+1(cos(φij) +13)2, (5.6)is calculated for each molecule, and q distributions are obtained by averagingover ∼ 1000 configurations. In Equation (5.6), φij is the angle between linesconnecting the oxygen atom of a given molecule to the oxygen atoms of itsi-th and j-th nearest neighbours.In liquid water, there are local density fluctuations associated with fluc-tuations in tetrahedral order, and these are responsible for water’s character-istic anomalies [125]. Here we explore whether or not electric fields influencedensity and structural fluctuations. Following earlier work [125], we classifyevery water molecule as either high-q (H) or low-q (L) depending on whetherits q-value was larger or smaller than the configurational median q-value. Wethen calculate partial radial distribution functions, gαβ(r) = hαβ(r)+1, andstructure factors,Sαβ(k) = δαβ + 4piρ√xαxβ∫ ∞0drr sin(kr)k hαβ(r), (5.7)where xα is the mole fraction of “species” α. At k = 0, structure factors mea-sure fluctuations: the density-density structure factor, SNN (k) = S(k), mea-sures density fluctuations, the concentration-concentration structure factor,SCC(k) = xHxL [xHSHH(k) + xLSLL(k)− 2√xHxLSHL(k)] , (5.8)measures concentration fluctuations, and, the density-concentration struc-995.3. Resultsture factor,SNC(k) = xHxL[SHH(k)− SLL(k) + 2xHxL√xHxLSHL(k)]. (5.9)measures the coupling between density and concentration fluctuations [126].The presence of ice-like molecules is determined using the CHILL algo-rithm [75]. The CHILL algorithm distinguishes between cubic, hexagonal,and interfacial ice-like molecules. An ice-like molecule is defined to belongto an ice cluster if it is within 0.33 nm [approximate location of the firstminimum in the oxygen-oxygen g(r)] of another ice-like molecule.5.3 Results5.3.1 Melting TemperaturePrevious simulations have demonstrated that upon supercooling water mod-els freeze quickly when interacting with a uniform electric field on the orderof 1 V/nm [47, 58, 62]. Consistent with those results, water interacting witha field of 1 V/nm froze spontaneously within 40 ns in all ten independentsimulations we ran at 270 K with 8000 molecules. According to classicalnucleation theory, ice is more likely to nucleate in deeply supercooled waterbecause the critical nucleus is smaller. If electric fields increase water’smelting temperature (TM ), applying a field would increase the degree ofsupercooling, ∆TS = TM −T , at a given temperature, and explain the rapidfreezing observed.There are several methods commonly used to estimate melting temper-1005.3. Resultsatures from molecular dynamics simulations [72, 101, 127]. In the directcoexistence method [72, 101, 127], NPT simulations are performed withliquid and solid phases in contact and the melting temperature is deter-mined by finding the lowest temperature at which the solid phase melts. Inmost cases, the initial interface is constructed, rather than spontaneouslygenerated in a low-temperature simulation. Using the direct coexistencemethod, Abascal et al. [101] found melting temperatures of 287 and 289 Kfor the six-site water model, with 865 and 1536 molecules, respectively.Field magnitude TM (K) ∆Sfus < mx >(L) < mx >(S)(V/nm) N = 8000 (J/mol K) at TM at TM0.0 289(Ref. [101]) 24(Ref. [72]) 0.00 0.00270.8 307.8 ± 2.4 33 0.65 0.981.0 313.5 ± 2.4 33 0.71 0.981.5 325.2 ± 3.3 32 0.79 0.982.0 333.2 ± 2.9 32 0.83 0.98Table 5.1: Melting temperature estimates for different fields obtained fromsimulations with 8000 molecules. Estimates of ∆Sfus are obtained fromNPT simulations with 1000 molecules. The last two columns are the averagepolarization for liquid and solid systems. L and S stand for liquid and solidrespectively.We use the direct coexistence method to estimate the melting tempera-ture of water interacting with an electric field. However, instead of placingliquid and constructed solid phases in contact, we simply take “half-frozen”configurations from simulations with a field of 1 V/nm at 270 K, as initialstates for the direct coexistence method. In the half-frozen configurationsused, between 40 and 50 % of the molecules had ice-like structure, as1015.3. ResultsFigure 5.1: Example starting configuration containing ∼ 4000 ice-likemolecules. The oxygen atoms of cubic ice-like molecules are shown in darkblue, hexagonal ice-like molecules in light blue, and liquid molecules in red.All hydrogen atoms are white.1025.3. Resultsidentified by the CHILL algorithm [75]. In comparison with configurationsgenerated by placing solid and liquid phases in contact, we would expectconfigurations from freezing trajectories to have lower interfacial energy, andan ice surface (or surfaces) well suited for crystal growth. A sample startingconfiguration is shown in Fig. 5.1. The melting temperatures, averagedover four independent simulations with N = 8000, are reported in Table5.3.1 for different field magnitudes. The melting temperature increases withfield magnitude. For example, the melting temperature increases by ∼ 24K with a field of 1 V/nm, and by ∼ 44 K with a field of 2 V/nm, assuming289 K is the melting temperature in the absence of an electric field [101].Melting temperatures determined using smaller simulations (N = 1000) areconsistent with those reported for larger systems and show a similar trend.We note that to avoid possibly biasing the crystal structure in freezingsimulations, it is desirable to allow the box lengths to fluctuate indepen-dently, as was done in earlier simulations involving hexagonal ice [128].However, in the present simulations with an external field, the box geometryfluctuated wildly if the box lengths were allowed to fluctuate independently.Therefore, in our simulations the box length ratios were held fixed. Thisis not an issue for cubic ice, which is the stable form of fully polarized ice.Moreover, we use a relatively large simulation cell (N = 8000), and initialconfigurations taken from freezing trajectories, which should mitigate anyinfluence of the imposed cell constraints.With a field of 1 V/nm, the six-site water model freezes readily at 270K (∆TS = 44 K). At the same field, ten independent 50 ns simulationswith 8000 molecules showed no sign of freezing at 280 K (∆TS = 34 K).1035.3. ResultsThus, model water interacting with an electric field does not freeze on thetimescales of our simulations until ∆TS ∼ 40 K, similar to the supercoolingrequired to freeze real water on experimental timescales. Our results suggestthat electric fields induce freezing in simulations of supercooled water mainlybecause the true melting point is increased by the field. This means thatthe true ∆TS is considerably larger than one would expect based on themelting point of water in the absence of a field. However, the nucleationbarrier remains high, and the liquid must be supercooled by ∼ 40 K relativeto the true field-dependent melting point in order to freeze on simulationtimescales.5.3.2 Ice Clusters and NucleationIn this section we compare nucleation in polarized water with what is knownabout nucleation in unpolarized water. In addition, we examine ice clusterformation and ice growth in systems interacting with a field of 1 V/nm, aswell as in systems without an external field.Regardless of field, all supercooled liquid systems studied contained ice-like molecules with both cubic and hexagonal structure. However, nearlyall molecules in larger (> 10 molecules) ice clusters had cubic structurein presence of an electric field, with hexagonal structure appearing onlynear the cluster surface. Both cubic and hexagonal ice consist of repeatinghexagonal rings: hexagonal ice contains boat and chair ring configurations,while cubic ice contains only chair configurations. Water dipoles cannotperfectly align (i.e., < mx >= 1) in the boat configuration, but can in thechair configuration; thus, only cubic ice grows in strongly polarized water.1045.3. ResultsNumber of Average number LargestT ∆TS Field independent Time of ice-like cluster(K) (K) (V/nm) simulations (ns) molecules observed260 29 0 10 40 33 11270 19 0 1 40 17 5270 28 0.8 1 40 24 14270 44 1.0 10 40 43 45280 34 1.0 10 50 17 11290 24 1.0 9 50 11 8Table 5.2: Summary of ice-like molecules and ice clusters observed indifferent systems; configurations from each simulation were saved at 10ps intervals. For some systems, multiple independent simulations wereperformed as indicated. At 270 K, freezing occurred within 40 ns for eachof the ten simulations performed with a field of 1 V/nm; the tabulated datawere taken from pre-freezing configurations. For all systems, the largestcluster reported is the largest ice cluster (see text for definition of cluster)observed that subsequently melted.We report the average numbers of ice-like molecules observed with dif-ferent field magnitudes and temperatures in Table 5.2, and show ice clusterdistributions for polarized (1 V/nm, 270 K and 280 K) and unpolarized (260K) systems in Fig. 5.2. All results are for simulations with 8000 molecules.The number of ice-like molecules increases with decreasing temperature, asdo the number and size of ice clusters. However, less than 1 % of watermolecules had ice-like structure in all liquid systems considered, even whensupercooled by more than 30 K. For example, 0.2 % of the molecules in asystem interacting with a 1 V/nm field have ice-like structure when super-cooled by 34 K (T = 280 K). In the absence of a field, liquid supercooled by29 K (T = 260 K) contains approximately twice as many ice-like molecules(0.4 %). Fewer ice-like molecules may form when interacting with a field at1055.3. Resultsa given ∆TS because dipole fluctuations are restricted in ice-like structures.Nevertheless, larger ice-like clusters, consisting of eight or more molecules,occurred with the same frequency in polarized and unpolarized systems, asshown in the inset of Fig. 5.2.-5 0 5 10 15 20 25 30 35 1  2  3  4  5  6  7  8Number of clusters per config.Number of molecules in cluster 0 0.001 0.002 6  8  10  12Figure 5.2: Average number of clusters (per configuration) of a given sizeobserved in NPT (1 bar) simulations with 8000 molecules. Results areshown for: 0 V/nm, 260 K (red); 1 V/nm, 280 K (blue); 1 V/nm, 290K(black). Standard errors estimated from ten independent simulations of 40ns (260 K) or 50 ns (280 and 290 K) are also shown.Although all systems examined contain very few ice-like molecules, thepositions of ice-like molecules are strongly correlated, as shown by the pairdistribution functions plotted in Fig. 5.3. Correlations between ice-likemolecules grow with decreasing temperature, but are stronger in systemsinteracting with an electric field, even when they contain fewer ice-likemolecules. However, most ice-like molecules form in isolation, and clustersconsisting of more than two or three molecules are rare in all systems (Fig.1065.3. Results5.2).-0.5 0 0.5 1 1.5 2 2.5 3 3.5 4 0.2  0.4  0.6  0.8  1  1.2  1.4  1.6  1.8  2g(r)r (nm)NPT; 0.0 V/nm; 260K;NPT; 1.0 V/nm; 270K;NPT; 1.0 V/nm; 280K;Figure 5.3: Oxygen-oxygen radial distribution functions for ice-likemolecules obtained from simulations (N = 8000, P = 1 bar) for threesystems, as indicated in the legend.Ten independent simulations with 8000 molecules, and four with 32000molecules, were carried out at 270 K (∆TS ∼ 44 K) with a field of 1 V/nm. Inall cases freezing was observed within 40 ns. Before freezing, the largest icecluster that formed and then subsequently melted contained ∼ 45 molecules.Thus, larger clusters freeze, and we estimate that the critical nucleus at 270K contains ∼ 50 molecules.With a field of 1 V/nm, spontaneous freezing did not occur for anytemperature considered above 270 K. Therefore, to estimate the size ofcritical nuclei at higher temperatures, we used configurations from the large(N = 32000) freezing simulations at 270 K as starting configurations insimulations at higher temperatures. The selected configurations each con-1075.3. ResultsFigure 5.4: Example starting configuration with an ice “nucleus” of ∼ 500molecules. Only cubic (blue) and hexagonal (red) ice-like molecules areshown.1085.3. Resultstained one large ice cluster of between 480 and 550 molecules along withother smaller ice clusters and isolated ice-like molecules, as summarized inTable 5.3. An example starting configuration is shown in Fig. 5.4. Theice clusters grew at 290 K (∆TS ≈ 24 K), but melted at 295 K (∆TS ≈ 19K) for all four initial configurations. Thus, the critical nucleus shrinks from∼ 500 molecules when ∆TS ≈ 24 K to ∼ 50 molecules when ∆TS ≈ 44 K,in water interacting with a field of 1 V/nm.The size of the critical nucleus has not been previously reported at anytemperature for the six-site model, in the absence of an electric field. Todetermine whether electric fields affect the size of the critical nucleus for agiven ∆TS , we used the same four starting configurations described above(Table 5.3), each containing an ice cluster of ∼ 500 molecules, in simulationswithout an external electric field. The ice clusters melted within 5 ns at 265K (∆TS ≈ 24 K), but unpolarized ice grew on the cluster surfaces at 260 K(∆TS ≈ 29 K), in all simulations. Growth of unpolarized ice on a polarizednucleus has been observed in previous simulations [115, 116] (Chapter 3 and4). While we have not precisely determined the lowest temperature at whicha 500-molecule ice cluster grows, growth occured with a 1 V/nm field when∆TS ≈ 24 K and without a field at ∆TS ≈ 29 K. The two temperaturesare reasonably close, suggesting that the size of the critical nucleus for thesix-site model depends mostly on the degree of supercooling. This furthersupports the conclusion noted above, that the effect of the field comes mainlyfrom its influence on the melting point.A critical nucleus of ∼ 500 molecules when ∆TS ≈ 29 K (no field case)is somewhat smaller than estimates for the TIP4P/2005 and TIP4P/Ice1095.3. ResultsNumber of clusters Size of Size ofSimulation with >10 largest 2nd largest NC NH NImolecules cluster cluster1 2 485 11 583 110 11832 2 515 12 619 102 11713 2 507 11 593 127 11734 2 545 31 646 111 1232Table 5.3: Summary of ice-like molecules in starting configurations used toestimate the critical ice nucleus size. The starting configurations were takenfrom simulations performed with 32000 water molecules interacting with afield of 1 V/nm at 270 K, where freezing eventually occurred. The numberof molecules meeting the condition for interfacial ice (NI), as identified bythe CHILL algorithm, are recorded in the table, but only cubic (NC) andhexagonal (NH) ice-like molecules were included in determining the size ofice clusters. These four starting configurations were used in simulations toestimate the highest temperature at which the initial ice cluster nucleatesice, both with and without a field.models (∼ 600 molecules when ∆TS = 35 K) [120]. While the discrepancymight well be due to model differences, our starting configurations were verydifferent from those used by Sanz et al. [120], which contained sphericalclusters cut from a hexagonal crystal. In our simulations, the startingconfigurations 1) contained ice clusters that had cubic rather than hexagonalstructure, 2) were taken from freezing trajectories, and 3) were polarized,denser, and contained other ice-like molecules, in addition to the ∼ 500molecule ice cluster. With respect to the last point, ice growth was relativelyslow at 260 K, taking ∼ 5-10 ns for the cluster to grow to 1000 molecules,allowing sufficient time for system relaxation. For example, in one simulationwith 32000 molecules, an initial cluster containing 485 molecules shrank to448 molecules, and the total number of ice-like molecules decreased by nearly1105.3. Results100 molecules, before irreversible ice growth occurred. The fact that the icecluster decreased in size suggests that the critical ice nucleus is somewhatsmaller than our estimate.5.3.3 Thermodynamic FunctionsFig. 5.5 shows how the average polarization, configurational energy, field-interaction energy, entropy [Equation (5.4)], and the Gibbs free energychange with field strength. Water becomes more polarized as the fieldincreases; as a result, water molecules interact more favourably with thefield, and lose rotational entropy. The average configurational energy alsodecreases as water molecules align with the field, indicating the presenceof more energetically favourable structures. The favourable configurationaland field interaction energies compensate for the loss of entropy, and thefree energy of liquid water decreases with increasing field (Fig. 5.5).Usually, liquids are stable with respect to solids because of their greaterentropy. Liquid water loses entropy when interacting with a uniform electricfield, offering a possible explanation for the increased melting temperature.However, liquid water becomes more stable when interacting with a field(∆G < 0), and therefore the melting temperature can only increase if icealso becomes more stable. Table 5.3.1 shows that in ice the molecular dipolesalign essentially perfectly with the field, whereas significant orientationaldisorder persists in liquid water at the freezing transition. The entropyof fusion is also larger in the presence of an electric field, indicating thatice loses slightly more entropy than liquid water, likely due to the higherpolarization. Although associated with a greater loss of rotational entropy,1115.3. Results0 1 2 3 4-15-10-50Energy (KJ/mol)Configurational energy differenceWater-field interaction energyGibbs Free energy differenceT∆ S0 1 2 3 4Electric field (V/nm)00.40.8Dipole order parameterFigure 5.5: Changes in thermodynamic properties (top) and the dipole orderparameter (bottom) as functions of field strength.1125.3. Resultsnear perfectly aligned dipoles give ice an energetic advantage through thefield-interaction energy. It is in fact this energetic advantage that increasesthe stability of ice with respect to liquid water, resulting in an increasedmelting temperature.To be more explicit, liquid water loses less entropy than ice due to fieldinteractions, and its configurational energy decreases with increasing field(Fig. 5.5), whereas the configurational energy of ice is nearly independentof field. Yet, liquid water loses stability with respect to ice because thepartially aligned dipoles in the liquid interact less favourably with the fieldthan the fully aligned dipoles in cubic ice. This suggests that in liquid waterthe molecular dipoles cannot perfectly align (and take full advantage ofthe favourable field-interaction energy) without losing entropy due to otherstructural changes.It is of interest to compare the observed entropy loss with the expectedideal entropy loss for a non-interacting dipolar fluid polarized by an electricfield. The partition function QID(E) for ideal, non-interacting, dipolarmolecules in the presence of a field is [129]QID(E) = QID(0)sinh yy , (5.10)where y = µwE/kT , E = |E|, µw is the permanent dipole moment for thesix-site water model, k is the Boltzmann constant, and QID(0) is the parti-tion function in the absence of a field. The change in the ideal Helmholtz1135.3. Resultsfree energy due to interaction with an applied field is given by∆AID(E) = AIDfield −AIDno field = −kT ln[QID(E)]. (5.11)The average polarization, P , for an ideal system at a specific temperatureis calculated through the Langevin function [129]P = L(y) = coth y − 1y . (5.12)The ideal entropy change ∆SID(E) can then be obtained∆SID(E) = SIDfield − SIDno field =∆U ID(E)−∆AID(E)T , (5.13)where U ID(E) = −µwEP , and AID(E) = −kT ln[QID(E)]. For our simu-lated NPT systems, we determined the change in entropy due to polariza-tion using Equation (5.4).Fig. 5.6 shows the entropy change as a function of polarization for idealand simulated systems at 300 K. The simulated system loses slightly lessentropy than the ideal dipolar liquid with the same average polarization,likely due to correlations present in the absence of an electric field. However,as the polarization approaches 1, (& 0.87), liquid water at 300 K loses moreentropy than the available ideal rotational entropy, indicating that liquidwater adopts a less entropically favoured structure when highly polarized.1145.3. Results0 0.2 0.4 0.6 0.8 1Polarization-20-100Entropy change (J/mol)NPT-systemtheoretical-predictionFigure 5.6: Entropy change as a function of polarization (< mx >) at 300K. The red circles are from Equation (5.13), and the blue triangles denotethe simulation results. The lines are to guide the eye.1155.3. Results5.3.4 Structural Changes with FieldOne might have expected polarization to significantly alter the structure ofliquid water. However, consistent with the largely ideal entropy changes dis-cussed above, the structure of liquid water remains qualitatively unchangedby polarization if < mx >. 0.87. This is demonstrated by tetrahedral orderparameter distributions (Fig. 5.7) and r2h(r) plots (Fig. 5.8). We note thatwater retains its local tetrahedral order (Fig. 5.7), and that hHH(r), hLL(r),and hHL(r) change very little when the liquid is polarized (Fig. 5.8). 0 0.001 0.002 0.003 0  0.2  0.4  0.6  0.8  1P(q)q270K; Ex=1.0 V/nm260K; Ex=0.0 V/nm270K; Ex=0.0 V/nmFigure 5.7: Comparison of tetrahedral order parameter distributions atconstant pressure for different temperatures and external fields, as indicatedin the legend.Despite these striking similarities, electric fields do induce small quanti-tative changes in the structure of water. In particular, we note that the localstructure becomes more tetrahedral (Fig. 5.7), and the radial correlationsincrease slightly (Fig. 5.8) in presence of a field. The structure of waterinteracting with an electric field is similar to that of slightly colder (by ∼ 51165.3. Results-0.05 0 0.05 0.1 0.15 0.2 0.25r2 h LL(r)270 K; 0 V/nm260 K; 0 V/nm280 K; 1.0 V/nm270 K; 1.0 V/nm-0.1-0.05 0 0.05 0.1 0.15 0.2 0.25 0.3r2 h HH(r)-0.05 0 0.05 0.1 0.15 0.2r2 h HL(r)-0.1-0.05 0 0.05 0.1 0.15 0.2 0.251 2r2 h(r)r (nm)Figure 5.8: Total (bottom) and partial (as labeled) pair correlation functionsat the temperatures and fields indicated in the legend. Multiplication by r2magnifies the structural features.1175.3. ResultsK) water in absence of a field, and is consistent with a more negative con-figurational energy (Fig. 5.5). However, the density of water increases withpolarization, for example, from 0.992 to 1.002 g/cm3 at 270 K (Table 5.4).Increased density is usually associated with less tetrahedral structure andwarmer water, thus, polarization allows water to accommodate energeticallyfavoured tetrahedral structures in a denser liquid. Nevertheless, the densityof polarized water decreases with decreasing temperature, consistent withwater’s characteristic density anomaly.T (K) Field (V/nm) Density (g/cm3)260 0 0.973270 0 0.992270 1.0 1.002280 1.0 1.019290 1.0 1.023Table 5.4: Densities in polarized and unpolarized water from simulationswith 8000 molecules at 1 bar.At higher fields (< mx > & 0.87 at 300 K) significant structural changesdo occur, with water losing some tetrahedral order and gaining long-rangecorrelations. Although loss of tetrahedral structure is usually accompaniedby a less favourable configurational energy, highly polarized water has alower configurational energy. The increased long-range correlations makethis structure less entropically favourable, perhaps accounting for the addi-tional loss of entropy, beyond that expected from the loss of rotational en-tropy that occurs when water is highly polarized. Water froze spontaneouslyin some NPT simulations at 300 K with a field of 5 V/nm; however, a high1185.4. Summary and Conclusionsdensity ice, rather than cubic ice, was formed. Nevertheless, cubic ice doesnot melt at such high fields and has both a lower configurational energyand a lower field-interaction energy than the high density form, suggestingthat water freezes to the high density ice structure because it is kineticallyaccessible, rather than thermodynamically stable.The concentration-concentration, SCC(k), and number-concentration,SNC(k), structure factors can be calculated from the partial structure factorsfor H and L molecules, as discussed in Section 5.2. In order to accuratelydetermine the structure factors at low wavenumbers, the relevant radialdistribution functions are corrected for finite-size effects following Salacuseet al. [130]. The magnitudes of S(0), SCC(0) and SNC(0) increase in thepresence of an electric field, indicating that both density and concentrationfluctuations increase, as does the coupling between them. The values ofS(0), SCC(0) and SNC(0) at 270 K with a field of 1.0 V/nm are comparableto those measured at 260 K in an unpolarized system, as shown in Fig. 5.9.Thus, apart from its higher density, polarized water appears colder, withmore tetrahedral local order, stronger structural correlations, and largercoupled density and concentration fluctuations.5.4 Summary and ConclusionsIn this chapter we show that uniform fields on the order of 1-2 V/nmsignificantly increase the melting point of the six-site water model. Theincreased stability of polarized cubic ice can be traced to its near perfectpolarization, whereas the coexisting liquid achieves only ∼ 80% polarization.1195.4. Summary and Conclusions-0.4-0.2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.80 10 20 30 40 50 60 70 80S CC(k); S(k); SNC(k)k (nm-1)SCC(k)SNC(k)S(k)270 K; 0.0 V/nm260 K; 0.0 V/nm280 K; 1.0 V/nm270 K; 1.0 V/nmFigure 5.9: Structure factors for the temperatures and electric fieldsindicated in the legend.1205.4. Summary and ConclusionsThus, the average field-dipole interaction energy favours the ice. Therefore,compared to the zero field case, at a given temperature polarized water ismore deeply supercooled relative to its true field-dependent melting point.This reduces the size of the critical nucleus, and helps explain why polarizedwater freezes readily on simulation timescales. We find spontaneous freezingat ∆TS = TM − T ≈ 40 K. We also show that the size of the critical icenucleus is mostly determined by ∆TS , regardless of field.Our observations also offer an explanation for heterogeneous ice nucle-ation by local electric fields. As noted above, uniform electric fields aslarge as those considered in this work would cause dielectric breakdown inreal water. However, water can remain stable when interacting with largelocal fields generated by ions or surface charge groups [55–57, 63]. Fieldsgenerated by some ice-nucleating materials may cause nearby water dipolesto partially align with the field, thereby effectively increasing ∆TS locally.Since the size of the critical nucleus decreases with increasing ∆TS , moleculesaffected by the local field can freeze, and nucleate bulk ice, as demonstratedin Chapter 3 and 4. Weaker electric fields may also nucleate ice, but at arate that is too slow to observe in simulations currently feasible.Although supercooled water contains rather few ice-like molecules (<1%), ice-like molecules associate strongly with each other, both with andwithout an interacting electric field. For the same degree of supercooling(∆TS), water molecules form ice-like structures less readily in presence ofan electric field, but the ice-like molecules that do form associate with eachother more strongly.It is interesting to note that the structure of water in a uniform electric1215.4. Summary and Conclusionsfield shows little qualitative change until very high polarization (& 87%)is achieved. At lower polarizations water maintains its normal tetrahedral,hydrogen-bonded structure and characteristic spatial correlations. Even thedensity anomaly and the anomalous increase in the structure factor at lowwavenumber are left intact by the field. Thus, atomic-level informationobtained from simulations with electric fields may provide some insight intohomogeneous ice nucleation in absence of a field.122Chapter 6Conclusions and Perspective6.1 Simulation Results on Ice Nucleation andElectrofreezingThis thesis systematically investigates the process of ice nucleation inducedby electric fields, and why external electric fields promote the freezing ofliquid water. The microscopic mechanism and the growth pattern of icenucleation are focuses of the investigation. Molecular dynamics simulationsare employed to simulate the occurrence of ice, and calculate structuraland dynamical properties, such as density and dipole order profiles, radialdistribution functions, ice-like molecule detections, diffusion coefficients, andso forth.In Chapter 2, models of water and the confinement geometries are de-scribed in detail, as is the model of electric fields used in this thesis. Thecomputational approaches of molecular dynamics simulation are introduced,and methods of data analysis of the simulation results are discussed.In Chapter 3, ice formation induced by a decaying field, which actsnear the entire surface, is described. It is demonstrated that a surface fieldcan create an effective ice nucleus in models of supercooled liquid water.1236.1. Simulation Results on Ice Nucleation and ElectrofreezingThe field must polarize only a very thin water layer (∼ 10 A˚), and thefield strength required is realistic on the relevant length scale. Ice withferroelectric cubic structure nucleates near the surface, and dipole disorderedcubic ice grows outward from the surface. Moreover, two different watermodels (the six-site and TIP4P/Ice models) are compared, and in both casesit is shown that a surface field can serve as a very effective ice nucleationcatalyst in supercooled water. We examine the influence of temperatureand two important field parameters, the field strength and distance fromthe surface over which it acts, on the nucleation process. For the six-sitemodel the highest temperature where we observe field-induced ice nucleationis 280 K (normal melting point is 289 K), and for TIP4P/Ice 270 K (normalmelting point is 270 K). The minimum electric field strength required tonucleate ice on the timescales of our simulations depends a little on how farthe field extends from the surface. If it extends 20 A˚, then a field strength of1.5 × 109 V/m is effective for both models. If the field extent is 10 A˚, thenstronger fields are required (2.5 × 109 V/m for TIP4P/Ice, and 3.5 × 109V/m for the six-site model). Our results demonstrate that fields of realisticstrength, that act only over a narrow surface region, can effectively nucleateice at temperatures not far below the freezing point. This supports thepossibility that local electric fields can be a significant factor influencingheterogeneous ice nucleation in physical situations.In Chapter 4, the surface field is modified so that it decays in twodirections. Field “bands” of different geometry all nucleate ice, providedthat the band is sufficiently large. Rectangular bands are very efficientif the width and thickness are & 0.35 nm, and & 0.15 nm, respectively,1246.1. Simulation Results on Ice Nucleation and Electrofreezingand the necessary dimensions are comparable for other geometries. Carefulanalysis of different systems reveals that ice strongly prefers to grow at(111) planes of cubic ice. This agrees with an earlier theoretical deductionbased on considerations of water-ice interfacial energies [111]. It is foundthat ice nucleated by field bands usually grows as a mixture of cubic andhexagonal ice, consistent with other simulations of ice growth [73, 75, 99,119], and with experiment [103, 106]. The mixture of ice Ih and Ic contrastswith simulations carried out with nucleating fields that span the entiresurface area, where cubic ice dominates, and hexagonal layers are very rarelyobserved [115, 116]. This discrepancy is mainly a simulation artifact relatedto finite sample size and periodic boundary conditions.We have found that ice nucleates efficiently when induced by an electricfield. However, the question of how an electric field enhances ice nucleationis also of interests. In Chapter 5, we investigate the structural and ther-modynamic changes that promote freezing when liquid water interacts witha uniform, external field. The melting temperature increases significantlywhen electric fields are applied. Fields of 1 ×109 V/m and 2 ×109 V/mincrease the melting point by 24 K and 44 K, respectively. The increasedmelting point is mainly due to the favourable interaction of near perfectlypolarized cubic ice with the applied field. For a fixed temperature, wedemonstrate that the size of the critical ice nucleus decreases with fieldpresent in the system mostly because the melting point, and hence the truedegree of supercooling, is increasing with field. On simulation timescales,ice nucleation is observed at 40 K below the field-dependent melting point,independent of the particular value of the field applied. We find that even1256.2. Future Workquite highly polarized liquid water retains its characteristic local structure,and the related anomalous properties of water. These results are relevantto the mechanism of heterogeneous ice nucleation by local surface fields asdiscussed in earlier chapters. Local fields can effectively increase the degreeof supercooling of locally polarized liquid, decreasing the size of the criticalnucleus in the region influenced by the field, hence facilitating ice nucleation.6.2 Future WorkOur simulation results support the idea that local electric fields could playa major role in heterogeneous ice nucleation, particularly for rough particleswith many surface structure variations, such as kaolinite. These particlescan serve as ice nuclei in environmentally realistic situations. One area offuture work would be the investigation of ice nucleation on more realisticsurfaces with different arrangements of charge. Replacing the smooth wallswith an atomistic surface would be a step in that direction. An atomisticsurface could be “decorated” with different charges, which might have alarge influence on ice nucleation near the surface. It has been reported thatwater freezes differently on positively and negatively charged surfaces inexperiments [57], and using a more realistic atomistic surface is one way tosimulate similar physical environments.Simulations of ice nucleation on realistic particle models are of greatinterest. Silver iodide is one of the best IN known [7] and can nucleate ice attemperatures as warm as −3 ◦C. Under atmospherically relevant conditions,two phases of silver iodide exist [131], β-AgI and γ-AgI, with either silver1266.2. Future Workor iodide ions exposed to the adjacent medium. The AgI faces are bothmade of chair conformation hexagonal rings of alternating silver and iodideions forming a bilayer, such that either all silver or all iodide ions occupy thetopside layer. The interaction between water and surfaces may be influencedby the Coulombic attraction between water and the positive silver or thenegative iodide ions. It is speculated that the “bonding” between watermolecules and different AgI faces is important in ice nucleation. Recentprogress has been made in this area in our research group [132].Many other surfaces are of much interest as possible ice nuclei. In biologi-cal systems, some proteins have been observed to promote ice nucleation [55].Simulations with water and proteins at different temperatures is anotherpromising area of research in ice nucleation.127Bibliography[1] B. J. Murray, D. OSullivan, J. D. Atkinsona and M. E. Webbb, Icenucleation by particles immersed in supercooled cloud droplets. Chem.Soc. Rev. 41, 6519, (2012).[2] J. H. Seinfeld and S. N. Pandis, Atmospheric chemistry and physics -From air pollution to climate change. 1326 pp., (John Wiley & Sons,Inc., New York, 1998).[3] D. A. Hegg and M. B. Baker, Nucleation in the atmosphere. Rep.Prog. Phys. 72, 056801, (2009).[4] M. B. Baker, Cloud microphysics and climate. Science, 276, 1072,(1997).[5] M. B. Baker and T. Peter, Small-scale cloud processes and climate.Nature, 451, 299, (2008).[6] B. Zuberi, A. K. Bertram, C. A. Cassa, L. T. Molina, and M. J.Molina, Heterogeneous nucleation of ice in (NH4)2SO2−H2O particleswith mineral dust immersions. Geophys. Res. Let. 29, 142, (2002).[7] H. R. Pruppacher and J. D. Klett, Microphysics of Clouds andPrecipitation (2nd rev. ed.; Kluwer Academic: Dordrecht, 1997).128Bibliography[8] S. Hartmann, D. Niedermeier, J. Voigtlander, T. Clauss, R. A.Shaw, H. Wex, A. Kiselev and F. Stratmann, Homogeneous andheterogeneous ice nucleation at LACIS: operating principle andtheoretical studies. Atmos. Chem. Phys., 11, 1753, (2011).[9] T. Bartels-Rausch, V. Bergeron, J. H. E. Cartwright, R. Escribano, J.L. Finney, H. Grothe, P. J. Gutie´rrez, J. Haapala, W. F. Kuhs, J. B. C.Pettersson, S. D. Price, C. I. Sainz-Dı´az, D. J. Stokes, G. Strazzulla, E.S. Thomson, H. Trinks and N. Uras-Aytemiz, Ice structures, patterns,and processes: A view across the icefields. Rev. Mod. Phys. 84, 885,(2012).[10] W. Cantrell and A. Heymseld, Production of ice in tropospheric cloudsa review. B. Am. Meteorol. Soc. 86, 795, (2005).[11] C. Hoose, J. E. Kristja´nsson, J-P. Chen and A. Hazra, A classical-theory-based parameterization of heterogeneous ice nucleation bymineral dust, soot, and biological particles in a global climate model.J. Atmos. Sci.67, 2483, (2010).[12] D. J. Cziczo, P. J. DeMott, C. Brock, P. K. Hudson, B. Jesse, S. M.Kreidenweis, A. J. Prenni, J. Schreiner, D. S. Thomson and D. M.Murphy, A method for single particle mass spectrometry of ice nuclei.Aerosol Sci. Technol., 37, 460, (2003).[13] P. J. DeMott, D. J. Cziczo, A. J. Prenni, D. M. Murphy, S.M. Kreidenweis, D. S. Thomson, R. Borys and D. C. Rogers,129BibliographyMeasurements of the concentration and composition of nuclei for cirrusformation. Proc. Natl. Acad. Sci. U. S. A., 100, 14655, (2003).[14] K. A. Pratt, P. J. DeMott, J. R. French, Z. Wang, D. L. Westphal, A.J. Heymsfield, C. H. Twohy, A. J. Prenni and K. A. Prather, In situdetection of biological particles in cloud ice-crystals. Nat. Geosci., 2,397, (2009).[15] J. M. Prospero, P. Ginoux, O. Torres, S. E. Nicholson and T. E.Gill, Environmental characterization of global sources of atmosphericsoil dust identified with the NIMBUS 7 Total Ozone MappingSpectrometer (TOMS) absorbing aerosol product. Rev. Geophys., 40,1002, (2002).[16] R. A. Glaccum and J. M. Prospero, Saharan aerosols over the tropicalNorth Atlantic-Mineralogy. Marine Geology 37(3-4), 295, (1980).[17] D. J. Cziczo, D. M. Murphy, P. K. Hudson and D. S. Thomson, Singleparticle measurements of the chemical composition of cirrus ice residueduring CRYSTAL-FACE. J. Geophys. Res.-Atm. 109, 4517, (2004).[18] M. Dymarska, B. J. Murray, L. Sun, M. L. Eastwood, D. A. Knopfand A. J. Bertram, Deposition ice nucleation on soot at temperaturesrelevant for the lower troposphere. J. Geophys. Res.-Atm. 111, D4,(2006).[19] M. L. Eastwood, S. Cremel, C. Gehrke, E. Girard and A. K. Bertram,Ice nucleation on mineral dust particles: Onset conditions, nucleationrates and contact angles. J. Geophys. Res.-Atm. 113, D22203, (2008).130Bibliography[20] F. Zimmermann, M. Ebert, A. Worringen, L. Schtz and S. Weinbruch,Environmental scanning electron microscopy (ESEM) as a newtechnique to determine the ice nucleation capability of individualatmospheric aerosol particles. Atmos. Environ. 41, 8219, (2007).[21] S. J. Cox, Z. Raza, S. M. Kathmann, B. Slater and A. Michaelides,The Microscopic Features of Heterogeneous Ice Nucleation May AffectThe Macroscopic Morphology of Atmospheric Ice Crystals. FaradayDiscuss., 167, 389, (2013).[22] J.F. Bell, III, Water on Planets. Proc. Int. Astron. Union, 15, 29,(2009).[23] J. H. E. Cartwright, B. Escribano, and C.I. Sainz-Dı´az, The mesoscalemorphologies of ice films: Porous and biomorphic forms of ice underastrophysical conditions. Astrophys. J., 687, 1406, (2008).[24] F. Franks, Nucleation of ice and its management in ecosystems. Phil.Trans. R. Soc. Lond. A, 361, 557, (2003).[25] R. E. Lee, G. J. Warren and L. V. Gusta, Biological Ice Nucleationand its Applications (American Phytopathological Society, St. Paul,1995).[26] K. E. Zachariassen and E. Kristiansen, Ice nucleation and antinucle-ation in nature. Cryobiology, 41(4), 257, (2000).[27] Y. Ma, L. Zhong, H. Zhang and C. Xu, Effect of applied electric field131Bibliographyon the formation and structure of ice in biomaterials during freezing.Proc.of the 10rd ICSD, 789, (2010).[28] D. Gurian-Sherman and S. E. Lindow, Bacterial ice nucleation:significance and molecular basis. FASEB J., 7(14), 1338, (1993).[29] S. E. Lindow, The role of bacterial ice nucleation in frost injury toplants. Annu. Rev. Phytopathol. 21, 363,(1983).[30] L. Mishchenko, B. Hatton, V. Bahadur, J. A. Taylor, T. Krupenkinand J. Aizenberg, Design of Ice-free Nanostructured Surfaces Basedon Repulsion of Impacting Water Droplets. ACS Nano, 4(12), 7699,(2010).[31] A. Alizadeh, M. Yamada, R. Li, W. Shang, S. Otta, S. Zhong, L.Ge, A. Dhinojwala, K. R. Conway, V. Bahadur, A. J. Vinciquerra,B. Stephens and M. L. Blohm, Dynamics of Ice Nucleation on WaterRepellent Surfaces. Langmuir, 28(6), 3180, (2012).[32] J. M. Strong-Gunderson, R. E. Lee, Jr., and M. R. Lee, TopicalApplication of Ice-Nucleating-Active Bacteria Decreases Insect ColdTolerance. Appl Environ Microbiol., 58(9), 2711, (1992).[33] R. E. Feeney, T. S. Burcham and Y. Yeh, Antifreeze glycoproteinsfrom polar fish blood. Ann. Rev. Biophys. Biophys. Chem., 15, 59,(1986).[34] Z. Gezgin, T. C. Lee and Q. Huang, Engineering functional nanothin132Bibliographymultilayers on food packaging: ice-nucleating polyethylene films. JAgric Food Chem., 61(21), 5130, (2013).[35] J. Li, M. P. Izquierdo and T. C. Lee, Effects of ice-nucleation activebacteria on the freezing of some model food systems. Int. J. Food Sci.Technol., 32, 41, (2003).[36] J. C. Liao and K. C. Ng, Effect of ice nucleators on snow making andspray freezing. Ind. Eng. Chem. Res., 29(3), 361, (1990).[37] A. Michaelides and K. Morgenstern, Ice nanoclusters at hydrophobicmetal surfaces. Nature Mater., 6, 597, (2007).[38] M. Volmer and A. Weber,A. Keimbildung in u¨bersa¨ttigten Gebilden.Z. Phys. Chem. (Leipzig), 119, 227 (1926).[39] L. Farkas, Velocity of nucleation in supersaturated vapours. Z. Phys.Chem. (Leipzig), 125, 236, (1927).[40] R. Becker and W. Do¨ring, Kinetische Behandlung der Keimbildung inu¨bersa¨ttigten Da¨mpfen. Ann. Phys., 24, 719, (1935).[41] J. Frenkel, Statistical theory of condensation phenomena. J. Chem.Phys., 7, 200, (1939).[42] J. Curtius, Nucleation of atmospheric aerosol particles. C. R.Physique, 7, 1027, (2006).[43] C. Hoose and O. Mu¨hler, Heterogeneous ice nucleation on atmosphericaerosols: a review of results from laboratory experiments. Atmos.Chem. Phys., 12, 9817, (2012).133Bibliography[44] C. A. Stan, S. K. Y. Tang, K. J. M. Bishop and G. M. Whitesides,Externally applied electric fields up to 1.6 × 105 V/m do not affecthomogeneous nucleation of ice in supercooled water. J. Phys. Chem.B, 115, 1089, (2011).[45] P. W. Wilson; K. Osterday and A. D. J. Haymet, The effects of electricfield on ice nucleation may be masked by the inherent stochastic natureof nucleation. CryoLetters, 30, 96, (2009).[46] H. M. Jones and E. E. Kunhards, The influence of pressure andconductivity on the pulsed breakdown of water. IEEETrans. DEI-1,1016, (1994).[47] I. M. Svishchev and P. G. Kusalik, Crystallization of liquid water in amolecular dynamics simulation. Phys. Rev. Lett. 73, 975, (1994).[48] G. R. Edwards and L. F. Evans, Effect of Surface Charge on IceNucleation by Silver Iodide. Trans. Faraday Soc. 58, 1649, (1962).[49] G. A. Dawson and G. R. Cardell, Electrofreezing of supercooledwaterdrops. J. Geophys. Res. 78, 8864, (1973).[50] J. B. Doolittle and G. Vali, Heterogeneous freezing nucleation inelectric fields. J. Atmos. Sci. 32, 375, (1975).[51] T. G. Gabarashvili and N. V. Gliki, Origination of the ice phase insupercooled water under the influence of electrically charged crystals.Atmos. Oceanic Phys. 3, 324, (1967).134Bibliography[52] M. A. Abbas and J. Latham, The electrofreezing of supercooled waterdrops. J. Meteor. Soc. Japan 47, 65, (1969).[53] H. R. Pruppacher, Electrofreezing of supercooled water. Pure Appl.Geophys. 104, 623, (1973).[54] B. A. Tinsley and G. W. Deen, Apparent tropospheric response toMeV-GeV particle flux variations: A connection via electrofreezing ofsupercooled water in high-level clouds? J. Geophys. Res. 96, 22283,(1991).[55] M. Gavish, J. L. Wang, M. Eisenstein, M. Lahav and L. Leiserowitz,The role of crystal polarity in alpha-amino acid crystals for inducednucleation of ice. Science. 256, 815, (1992).[56] I. Braslavsky and S. G. Lipson, Electrofreezing effect and nucleationof ice crystals in free growth experiments. Appl. Phys. Lett 72, 264,(1998).[57] D. Ehre, E. Lavert, M. Lahav and I. Lubomirsky, Water freezesdifferently on positively and negatively charged surfaces of pyroelectricmaterials. Science. 327, 672, (2010).[58] I. M. Svishchev and P. G. Kusalik, Electrofreezing of liquid water: amicroscopic perspective. J. Am. Chem. Soc. 118, 649, (1996).[59] W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey andM. L. Klein, Comparison of simple potential functions for simulatingliquid water. J. Chem. Phys. 79, 926, (1983).135Bibliography[60] H. J. C. Berendsen, J. R. Grigera and T. P. Straatsma, The missingterm in effective pair potentials. J. Phys. Chem. 91, 6269, (1987).[61] X. Xia and M. L. Berkowitz, The electric field induced restructuringof water at platinum/water interface: molecular dynamics computersimulation. Phys. Rev. Lett. 74, 3193, (1995).[62] R. Zangi and A. E. Mark, Electrofreezing of confined water. J. Chem.Phys. 120, 7123, (2004).[63] T. Croteau, A. K. Bertram and G. N. Patey, Observations of high-density ferroelectric ordered water in kaolinite trenches using MonteCarlo simulations. J. Phys. Chem. A, 114, 8396, (2010).[64] T. Croteau, A. K. Bertram and G. N. Patey, Adsorption and Structureof Water on Kaolinite Surfaces: Possible Insight into Ice Nucleationfrom Grand Canonical Monte Carlo Calculations. J. Phys. Chem. A.112, 10708, (2008).[65] D. R. Nutt and A. J. Stone, Ice nucleation on a model hexagonalsurface. Langmuir. 20, 8715, (2004).[66] X. L. Hu and A. Michaelides, Ice formation on kaolinite: lattice matchor amphoterism? Surf. Sci. 601, 5378, (2007).[67] X. L. Hu and A. Michaelides, Water on the hydroxylated (011) surfaceof kaolinite: from monomer adsorption to a flat 2D wetting layer. Surf.Sci. 602, 960, (2008).136Bibliography[68] F. Lu¨o¨nd, O. Stetzer, A. Welti and U. Lohmann, Experimental studyon the ice nucleation ability of size-selected kaolinite particles inimmersion mode. L. Geophys. Res. 115, D14201-1-14, (2010).[69] X. L. Hu and A. Michaelides, The kaolinite (001) polar basal plane.Surf. Sci. 604, 111, (2010).[70] L. Vrbka and P. Jungwirth, Homogeneous Freezing of Water Starts inthe Subsurface. J. Phys. Chem. B. 110, 18126, (2006).[71] J. L. F. Abascal, E. Sanz, R. Garcia Fernander, and C. Vega,A potential model for the study of ices and amorphous water:TIP4P/Ice. J. Chem. Phys. 122, 234511, (2005).[72] H. Nada and J. P. J. M. van der Eerden, An intermolecular potentialmodel for the simulation of ice and water near the melting point: Asix-site model of H2O. J. Chem. Phys. 118, 7401, (2003).[73] M. A. Carignano, Formation of stacking faults during ice growth onhexagonal and cubic substrates. J. Phys. Chem. C 111, 501 (2007).[74] G. P. Johari, On the coexistence of cubic and hexagonal ice between160 and 240 K. Phil. Mag. B 78, 375, (1998).[75] E. B. Moore, E. de la Llave, K. Welke, D. A. Scherlis and V. Molinero,Freezing, melting and structure of ice in a hydrophilic nanopore. Phys.Chem. Chem. Phys. 12, 4124, (2010).[76] J. C. Shelley and G. N. Patey, Boundary condition effects in137Bibliographysimulations of water confined between planar walls. Mol. Phys. 88,385, (1996).[77] M. P. Allen and D. J. Tildesley, Computer Simulations of Liquids(Clarendon: Oxford, U.K., 1987).[78] H. A. Lorentz, Ueber die Anwendung des Satzes vom Virial in derkinetischen Theorie der Gase. Ann. Phys. 248, 127, (1881).[79] D. Berthelot, Sur le Me´lange des Gaz. C. R. Hebd. Se´anc. Acad. Sci.126, 1703, (1898).[80] A. Y. Toukmaji and J. A. Board, Ewald summation techniques inperspective: A survey. Computer Phys. Communication. 95, 73,(1995).[81] M. Deserno and C. Holm, How to mesh up Ewald sums. I. A theoreticaland numerical comparison of various particle mesh routines. J. Chem.Phys. 109(18), 7678, (1998).[82] D. J. Evans, On the representation of orientation space. Mol. Phys.34, 317, (1977).[83] D. J. Evans and S. Murad, Singularity free algorithm for moleculardynamics simulation of rigid polyatomics. Mol. Phys. 34, 327, (1977).[84] R. J. Sadus, Molecular Simulation of Fluids: theory, algorithms, andobject-orientation. (Elsevier, New York, 1999).138Bibliography[85] I. P. Omelyan, On the numerical integration of motion for rigidpolyatomics: The modified quaternion approach. Comput. in Phys.12, 97, (1998).[86] D. J. Evans and G. P. Morriss, Statistical Mechanics of NonequilibriumLiquids, 2nd ed.; ANU E Press: Canberra, 2007; Chap. 5.2.[87] B. Hess, C. Kutzner, D. van Der Spoel and E. Lindahl, GROMACS 4:algorithms for highly efficient, load-balanced, and scalable molecularsimulation. J. Chem. Theory Comput. 4, 435, (2008).[88] S. J. Plimpton, Fast parallel algorithms for short-range moleculardynamics. J. Comput. Phys. 117, 1, (1995).[89] U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee, and L.G. Pedersen, A smooth particle mesh Ewald method. J. Chem. Phys.103, 8577, (1995).[90] R. W. Hockney and J. W. Eastwood, Computer Simulation UsingParticles. (CRC Press, 2010).[91] H. Kamberaj, R. J. Low, and M. P. Neal, Time reversible andsymplectic integrators for molecular dynamics simulations of rigidmolecules. J. Chem. Phys. 122, 224114, (2005).[92] W. C. Swope, H. C. Andersen, P. H. Berens, and K. R. Wilson,A computer-simulation method for the calculation of equilibrium-constants for the formation of physical clusters of molecules: Applica-tion to small water clusters. J. Chem. Phys. 76, 637, (1982).139Bibliography[93] S. Nose´, A molecular dynamics method for simulations in the canonicalensemble. Mol. Phys. 52, 255, (1984).[94] W. G. Hoover, Canonical dynamics: equilibrium phase-space distri-bution. Phys. Rev. A, 31, 1695, (1985).[95] B. Hess, H. Bekker, H. J. C. Berendsen, and J. G. E. M. Fraaije,LINCS: A linear constraint solver for molecular simulations. J. Comp.Chem. 18, 1463, (1997).[96] R. W. Hockney, S. P. Goel and J. Eastwood, Quiet high resolutioncomputer models of a plasma. J. Comp. Phys. 14, 148, (1974).[97] M. Parrinello and A. Rahman, Polymorphic transitions in singlecrystals: A new molecular dynamics method. J. Appl. Phys. 52, 7182,(1981).[98] S. Nose´ and M. L. Klein, Constant pressure molecular dynamics formolecular systems. Mol. Phys. 50, 1055, (1983).[99] A. V. Brukhno, J. Anwar, R. Davidchack and R. Handel, Challengesin molecular simulation of homogeneous ice nucleation. J. Phys.:Condense. Matter. 20, 494243, (2008).[100] T. Koop and B. Zobrist, Parameterizations for ice nucleation inbiological and atmospheric systems. Phys. Chem. Chem. Phys.11,10839, (2009).[101] J. L. F. Abascal, R. G. Ferna´ndez, C. Vega and M. A. Carignano, The140Bibliographymelting temperature of the six site potential model of water. J. Chem.Phys. 125, 166101-1-2, (2006).[102] R. G. Ferna´ndez, J. L. F. Abascal and C. Vega, The melting point ofice Ih for common water models calculated from direct coexistence ofthe solid-liquid interface. J Chem. Phys. 124, 144506-1-11, (2006).[103] B. J. Murray, D. A. Knopf and A. K. Bertram, The formation of cubicice under conditions relevant to Earth’s atmosphere. Nature, 434, 202,(2005).[104] G. P. Johari, Water’s size-dependent freezing to cubic ice. J. Chem.Phys. 122, 194504-1-5, (2005).[105] A. V. Gubskaya and P. G. Kusalik, The total molecular dipole momentfor liquid water. J. Chem. Phys.117, 5290, (2002).[106] T. L. Malkin, B. J. Murray, A. V. Brukhno, J. Anwar and C. G.Salzmann, Structure of ice crystallized from supercooled water. Proc.Natl. Acad. Sci. 109, 1041, (2012).[107] E. B. Moore and V. Molinero, Is it cubic? Ice crystallization fromdeeply supercooled water. Phys. Chem. Chem. Phys. 13, 20008,(2011).[108] J. C. Johnston amd V. Molinero, Crystallization, melting, and struc-ture of water nanoparticles at atmospherically relevant temperatures.J. Am. Chem. Soc. 134, 6650, (2012).141Bibliography[109] L. Tianshu, D. Donadio, G. Russo and G. Galli, Homogeneous icenucleation from supercooled water. Phys. Chem. Chem. Phys. 13,19807, (2011).[110] D. Quigley and P. M. Rodger, Metadynamics simulations of icenucleation and growth. J. Chem. Phys. 128, 154518-1-7, (2008).[111] T. Takahashi, On the role of cubic structure in ice nucleation. J.Crystal Growth, 59, 441, (1982).[112] J. P. D. Abbatt, Interactions of atmospheric trace gases with icesurfaces: adsorption and reaction. Chem. Rev. 103, 4783 (2003).[113] E. B. Moore and V. Molinero, Structural transformation in super-cooled water controls the crystallization rate of ice. Nature 479, 506(2011).[114] A. Reinhardt, J. P. K. Doye, E. G. Noya, and C. Vega, Local orderparameters for use in driving homogeneous ice nucleation with all-atom models of water. J. Chem. Phys. 137, 194504 (2012).[115] J. Y. Yan and G. N. Patey, Heterogeneous ice nucleation induced byelectric fields. J. Phys. Chem. Lett. 2, 2555 (2011).[116] J. Y. Yan and G. N. Patey, Molecular Dynamics Simulations of IceNucleation by Electric Fields, J. Phys. Chem. A. 116(26), 7057,(2012).[117] S. Bauerecker, P. Ulbig, V. Buch, L. Vrbka, and P. Jungwirth,Monitoring ice nucleation in pure and salty water via high speed142Bibliographyimaging and computer simulations. J. Phys. Chem. C 112, 7631(2008).[118] E. Pluhar˘ova´, L. Vrbka, and P. Jungwirth, The effect of surfacepollution on homogeneous ice nucleation: A molecular dynamics study.J. Phys. Chem. C 114, 7831 (2010).[119] P. Pirzadeh and P. G. Kusalik, On understanding stacking faultformation in ice. J. Am. Chem. Soc. 133, 704 (2010).[120] E. Sanz, C. Vega, J. R. Espinosa, R. Caballero-Bernal, J. L. F.Abascal and C. Valeriani, Homogeneous ice nucleation at moderatesupercooling from molecular simulations. J. Am. Chem. Soc. 135,15008, (2013).[121] L. R. Maki, E. L. Galyan, M. M. Chang-Chien and D. R. Caldwell,Ice Nucleation Induced by Pseudomonas syringae. Appl Microbiol. 28,456, (1974).[122] H. A. Kolodziej, B. P. Jones and M. Davies, High field dielectricmeasurements in water. J. Chem. Soc. Faraday Trans. II. 71, 269,(1975).[123] M. Orlowska, M. Havet, and A. Le-Bail, Controlled ice nucleationunder high voltage DC electrostatic field conditions. Food Res. Int.42, 879, (2009).[124] J. R. Errington and P. G. Debenedetti, Relationship between143Bibliographystructural order and the anomalies of liquid water. Nature, 409, 318,(2001).[125] S. D. Overduin and G. N. Patey, Understanding the structure factorand isothermal compressibility of ambient water in terms of localstructure environments. J. Phys. Chem. B. 116, 12014, (2012).[126] D. T. Bowron, J. L. Finney, A. Hallbrucker, I. Kohl, T. Loerting, E.Mayer and A. K. Soper, The local and intermediate range structures ofthe five amorphous ices at 80 K and ambient pressure: a Faber-Zimanand Bhatia-Thornton analysis. J. Chem. Phys. 125, 194502, (2006).[127] A. J. C. Ladd and L. V. Woodcock, Interfacial and co-existenceproperties of the Lennard-Jones system at the triple point. Mol. Phys.36, 611, (1978).[128] C. Vega, M. Martin-Conde, and A. Patrykiejew, Absence of superheat-ing for ice Ih with a free surface: A new method of determining themelting point of different water models. Mol. Phys. 104, 3583 (2006).[129] T. L. Hill, An Introduction to Statistical Thermodynamics. (Addison-Wesley, London, 1962).[130] J. J. Salacuse, A. R. Denton, and P. A. Egelstaff, Finite-size effectsin molecular dynamics simulations: Static structure factor andcompressibility. I. Theoretical method. Phys. Rev. E 53, 2382 (1996).[131] N. Hull and D. A. Keen, Pressure-induced phase transitions in AgCl,AgBr and AgI. Phys. Rev. B, 59(2), 750, (1999).144Bibliography[132] S. Zielke and G. N. Patey, Research of ice nucleation on AgI faces isin progress.145


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



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"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items