UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Molecular dynamics simulations of heterogeneous ice nucleation by atmospheric aerosols Zielke, Stephen Alexander 2016

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

Item Metadata

Download

Media
24-ubc_2016_november_zielke_stephen.pdf [ 95.41MB ]
Metadata
JSON: 24-1.0314173.json
JSON-LD: 24-1.0314173-ld.json
RDF/XML (Pretty): 24-1.0314173-rdf.xml
RDF/JSON: 24-1.0314173-rdf.json
Turtle: 24-1.0314173-turtle.txt
N-Triples: 24-1.0314173-rdf-ntriples.txt
Original Record: 24-1.0314173-source.json
Full Text
24-1.0314173-fulltext.txt
Citation
24-1.0314173.ris

Full Text

Molecular Dynamics Simulations of HeterogeneousIce Nucleation by Atmospheric AerosolsbyStephen Alexander ZielkeB.Sc. Chemistry, conc. Nanoscience, University of Calgary, 2012A THESIS SUBMITTED IN PARTIAL FULFILLMENTOF THE REQUIREMENTS FOR THE DEGREE OFDoctor of PhilosophyinTHE FACULTY OF GRADUATE AND POSTDOCTORALSTUDIES(Chemistry)The University of British Columbia(Vancouver)September 2016c© Stephen Alexander Zielke, 2016AbstractWater droplets in the atmosphere do not freeze homogeneously until -38◦C. Freez-ing at warmer temperatures requires heterogeneous ice nuclei (IN). Despite theimportance of ice in the atmosphere, little is known about the microscopic mech-anisms of heterogeneous ice nucleation. This thesis employs molecular dynam-ics simulations to investigate ice nucleation by silver iodide, kaolinite, potassiumfeldspar, gibbsite, and a protein.Silver iodide is one of the best known ice nucleating agents. We examinedseven surfaces of silver iodide and observed ice nucleation on three surfaces. Thesurfaces that nucleated ice organized the first layer of water molecules into a con-figuration resembling ice, such as chair conformed hexagonal rings. Surfacesthat do not nucleate ice do not organize water into icelike configurations, suchas planar rings. Results suggest lattice mismatch is insufficient in predicting icenucleation, and a finer atomistic match is required.Finite silver iodide disks and plates were used to probe the relationship be-tween the size of a nucleating surface and maximum temperature of ice nucle-ation. Larger disks nucleated ice at warmer temperatures than smaller disks byiiforming larger initial cluster of ice which could reach the critical size easier thanhomogeneously formed clusters.Kaolinite is a common clay known to nucleate ice. Our simulations investi-gated both sides of the (001) surface and found both sides able to nucleate ice.The Al-surface was simulated with varying degrees of freedom of motion. Anoptimum amount of movement was required to nucleate ice as the surface needsto adapt into a configuration favorable to ice. Ice nucleated on the Si-surface viathe formation of a novel composite surface structure which facilitated bulk icenucleation.Potassium feldspar simulations explored three variations of the two primarycleavage planes. All surfaces failed to nucleate ice and density profiles suggestthat the surfaces are unlikely to nucleate ice. We succeeded in nucleating ice ongibbsite with prepared surface conformations compatible with ice.Biological IN, such as ice nucleation proteins, are among the most efficient IN.We attempted to simulate ice nucleation via a protein, but were unable to achieveice nucleation.iiiPrefaceA large portion of this thesis contains previously published material which is re-produced here with permission. The three papers which contain the material are:• Zielke, S.A.; Bertram, A.K.; Patey, G.N. A Molecular Mechanism of IceNucleation on Model AgI Surfaces. J. Phys. Chem. B 2014 119, 9049-9055. [Chapter 3]• Zielke, S.A.; Bertram, A.K.; Patey, G.N. Simulations of Ice Nucleation byModel AgI Disks and Plates. J. Phys. Chem. B 2016 120, 2291-2299.[Chapter 4]• Zielke, S.A.; Bertram, A.K.; Patey, G.N. Simulations of Ice Nucleation byKaolinite (001) with Rigid and Flexible Surfaces. J. Phys. Chem. B 2016120, 1726-1734. [Chapter 5]Chapters 3, 4, and 5 come directly from the relevant papers. Minor changes havebeen made to these papers to make the whole thesis more cohesive. Portions ofChapters 1 and 2 come from the introductions and methods sections of the papers.Appendices A, B, and C are the supporting information from the papers.ivThe author of this thesis carried out all of the simulations, analysis, creation offigures and tables, and co-wrote the published papers. Conception of the researchwas done collaboratively by the authors of the papers.Chapters begin with a short quotation from the novel ‘Cat’s Cradle’ by KurtVonnegut.[Vonnegut, K. Cat’s Cradle. Dell Publishing, New York, 1963.] Theplot of the novel involves a phase of ice called ice-nine, which has a melting pointof 45◦C. Ice-nine was created by altering the way the water molecules are ar-ranged in the ice, resulting in the high melting point. Part of the idea for this camefrom Kurt’s brother, Bernard, who discovered silver iodide as an ice nucleus.[29]Simulations with silver iodide form several chapters of this thesis, thus I found itreasonable to include these quotes. Selection of some quotes were intended to behumorous. Please don’t read too much into them.vTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xivList of Symbols and Abbreviations . . . . . . . . . . . . . . . . . . . . . xliDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xlvii1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Kinetics of Ice Nucleation . . . . . . . . . . . . . . . . . . . . . . 21.2 Atmospheric Aerosols as Ice Nuclei . . . . . . . . . . . . . . . . 51.3 Ice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81.4 Ice Nucleating Materials . . . . . . . . . . . . . . . . . . . . . . 9vi1.5 Previous Simulation Studies of Ice Nucleation . . . . . . . . . . . 121.6 Objectives and Outline . . . . . . . . . . . . . . . . . . . . . . . 142 General Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 162.1 Molecular Dynamics . . . . . . . . . . . . . . . . . . . . . . . . 172.2 Interactions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 182.3 Equations of Motion . . . . . . . . . . . . . . . . . . . . . . . . . 252.4 Force Fields and Models . . . . . . . . . . . . . . . . . . . . . . 292.5 Simulation Setup . . . . . . . . . . . . . . . . . . . . . . . . . . 322.5.1 Geometry of AgI, Kaolinite, Feldspar, and Gibbsite SlabSimulations . . . . . . . . . . . . . . . . . . . . . . . . . 322.6 Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 353 A Molecular Mechanism of Ice Nucleation on Model AgI Surfaces . 393.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 393.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 403.3 Polarization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 413.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 433.4.1 Why do some model AgI surfaces nucleate ice and othersnot? . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 473.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 624 Simulations of Ice Nucleation by Model AgI Disks and Plates . . . . 634.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63vii4.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 644.3 Polarization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 674.4 Rational . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 684.5 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . 704.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 905 Simulations of Ice Nucleation by Kaolinite (001) with Rigid andFlexible Surfaces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 925.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 925.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 935.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 955.3.1 Varying Flexibility of the Al-Surface . . . . . . . . . . . . 955.3.2 Non-crystallographic Rigid Al-surfaces . . . . . . . . . . 1085.3.3 Si-surface . . . . . . . . . . . . . . . . . . . . . . . . . . 1135.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1196 Ice Nucleation Abilities of Potassium Feldspar and Gibbsite . . . . . 1216.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1216.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1226.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1246.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1337 Attempted Ice Nucleation by a Protein . . . . . . . . . . . . . . . . . 1347.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134viii7.2 Strategy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1357.3 Methods and Results . . . . . . . . . . . . . . . . . . . . . . . . 1377.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1428 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1438.1 Summary of Thesis . . . . . . . . . . . . . . . . . . . . . . . . . 1438.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 147Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148A Appendix to Chapter 3 . . . . . . . . . . . . . . . . . . . . . . . . . . 164A.1 Results for the TIP4P/Ice model . . . . . . . . . . . . . . . . . . 164B Appendix to Chapter 4 . . . . . . . . . . . . . . . . . . . . . . . . . . 172B.1 Comparison of force fields . . . . . . . . . . . . . . . . . . . . . 172C Appendix to Chapter 5 . . . . . . . . . . . . . . . . . . . . . . . . . . 181ixList of TablesTable 2.1 Parameters for water models used in this thesis. . . . . . . . . . 30Table 2.2 Parameters for the AgI and CLAFF[122] forcefields used inthis thesis. AgI values[35] are specific to the water model used. 31Table 3.1 Summary of the surfaces considered. XY dimensions, the num-bers of water molecules, and bulk liquid densities of liquid wa-ter are given. All surfaces except β -AgI(101¯0) were simulatedwith both silver and iodide ions exposed, and the surfaces are5 nm apart measured from the closest ions. . . . . . . . . . . . 41Table 3.2 Summary of the AgI faces considered, the associated mismatcheswith Ih and Ic, and the polymorph(s) of ice formed . . . . . . . 44Table 4.1 Summary of sytems considered. System C2 contains 2 mirror-ing plates, and system Dr is a rectangular plate. X/Y/Z givesthe dimensions of the simulation cell. ∆Ts reports the super-cooling required to nucleate ice. Densities are for the six-sitemodel at 300 K. . . . . . . . . . . . . . . . . . . . . . . . . . 66xTable 5.1 For each system considered, the XY cell dimensions, the num-ber of water molecules, and the bulk water densities at 300 Kare given. The Z cell dimension is 8.42661 nm in all cases (sur-faces are 6 nm apart measured from the edges of the unit cells),except for the large Si-surface, which is 1 nm larger. . . . . . . 94Table 5.2 Lattice mismatch percentages from equation 5.1. Numbers insquare brackets give the values of Im and In ([Im : In]) used.Numbers in parentheses are the lattice parameters. The onlygood match is highlighted in bold. . . . . . . . . . . . . . . . . 107Table 5.3 Prepared rigid surfaces. Yes and No indicate if ice nucleated ornot. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112Table 6.1 Simulation parameters and system properties. For each sys-tem considered, the XY cell dimensions, the number of watermolecules, and the bulk water densities at 300 K are given. . . . 124Table 6.2 Scalings of the gibbsite (001) surfaces explored. For each scal-ing the lattice parameters are given, along with the match to thelattice constants of ice Ih. . . . . . . . . . . . . . . . . . . . . 125Table A.1 List of simulations performed using the six-site model of waterat 270 K (supercooled by 20 K) . . . . . . . . . . . . . . . . . 165Table A.2 List of simulations performed using the TIP4P/Ice water modelat 250 K (supercooled by 20 K) . . . . . . . . . . . . . . . . . 166xiTable A.3 List of simulations performed to determine the highest temper-ature at which ice will nucleate on the Ag exposed β -AgI(0001)face . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 167Table A.4 List of simulations performed using non-standard charges. Thesix-site simulation with charges of ±0.8 e on the β -AgI(0001)surface only froze after 473 ns. . . . . . . . . . . . . . . . . . 169Table B.1 A list of simulations performed with the six-site model. Sys-tems marked with a ‘P’ were conducted at 1 bar pressure and a‘0’ indicates the Ag and I atoms were neutral. . . . . . . . . . . 176Table B.2 Continuation of Table B.1. . . . . . . . . . . . . . . . . . . . . 177Table B.3 A list of simulations performed with the TIP4P/2005 model . . 178Table C.1 Rigid-kaolinite simulations. . . . . . . . . . . . . . . . . . . . 181Table C.2 Free-H simulations. . . . . . . . . . . . . . . . . . . . . . . . 182Table C.3 Free-OH simulations . . . . . . . . . . . . . . . . . . . . . . . 182Table C.4 Free-OH simulations started from free-H simulations, as de-scribed in Chapter 5. TIP4P/Ice simulations were at 230 K andsix-site simulations at 240 K. Divisions in table separate simu-lations started from different Free-H simulations. . . . . . . . . 183Table C.5 Ice-prepared surface simulations. Simulations marked with anasterisk used ice-prepared surfaces produced by the other modelof water. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 184Table C.6 Liquid-prepared surface simulations. . . . . . . . . . . . . . . 184xiiTable C.7 “Ideal” surface simulations. . . . . . . . . . . . . . . . . . . . 185Table C.8 Vacuum-prepared surface simulations. . . . . . . . . . . . . . 185Table C.9 Si-surface simulations. The simulation marked with an asteriskis the large Si-surface simulation. . . . . . . . . . . . . . . . . 185xiiiList of FiguresFigure 1.1 Plot of eq 1.1 (black line). The red line is the surface energyterm and the blue line is the absolute value of the chemicalpotential term. . . . . . . . . . . . . . . . . . . . . . . . . . . 3Figure 2.1 Six-site and TIP4P/Ice water models. Oxygen atoms, hydro-gen atoms, lone pairs, and m-site are red, black, teal, andpink respectively. TIP4P/2005 shares the same basic shapeas TIP4P/Ice. . . . . . . . . . . . . . . . . . . . . . . . . . . 20Figure 2.2 A diagram of the simulation cells used in much of this the-sis. Opposing slabs of the material under consideration areshown in red, with the water between the slabs shown in blue.Vacuum areas at the end of the cell are white. . . . . . . . . . 34xivFigure 3.1 Polarization profiles for the six-site model. Panels A, B, andC show results for Ag exposed β -AgI(0001), Ag exposed γ-AgI(111), and Ag exposed γ-AgI(001), respectively. The red,green, and blue curves are the X, Y and Z components, re-spectively. Note the shifted scale in Panels A and B. . . . . . . 42Figure 3.2 Polarization profiles for the TIP4P/Ice model. Details are asin Fig. 3.1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43Figure 3.3 The number of ice molecules as a function of time detected bythe CHILL algorithm in simulations on three AgI surfaces: Agexposed β -AgI(0001) (red), Ag exposed γ-AgI(111) (green),and Ag exposed γ-AgI(001) (blue). The solid lines and dottedlines represent the number of Ic and Ih molecules, respectively. 45Figure 3.4 Snapshots of ice growth at different times. Panels A-C are forthe silver exposed β -AgI(0001) face at 0, 11 and 36 ns, re-spectively. Panels D-F are for the silver exposed γ-AgI(001)face at 0, 3 and 6 ns, respectively. Only the first layer ofeach AgI slab is shown. Silver and iodide ions are silver andgreen, respectively. Oxygen atoms of liquid water, Ic, and Ihare red, yellow and blue, respectively, and all hydrogen atomsare black. Note that panels D-F have been rotated to show the(110) plane of Ic, and this is why the water density can appearnonuniform in those snapshots. . . . . . . . . . . . . . . . . . 46xvFigure 3.5 Snapshots of ice growing on the Ag exposed γ-AgI(111) sur-face. Panels A-C are taken at 0 ns, 7 ns, and 15 ns, respec-tively. Silver and iodide ions are silver and green, respectively.Oxygen atoms of liquid water, Ic, and Ih are red, yellow andblue, respectively, and all hydrogen atoms are black. The pan-els are rotated as in Fig. 1D-E. . . . . . . . . . . . . . . . . . 47Figure 3.6 Snapshots of liquid water on β -AgI(0001) surfaces. Panel A isa top view of water (supercooled by 4 K) on the silver exposedface, and Panel C is a side view of a single hexagonal ring onthis surface. Panels B and D are corresponding views of water(supercooled by 20 K) on the iodide exposed face. Panel E is aside view of a chair structure from ice. Silver ions, iodide ions,oxygen atoms, and hydrogen atoms are silver, green, red andblack, respectively. The yellow lines highlight the hexagonalarrangement of the water on the surface. . . . . . . . . . . . . 49xviFigure 3.7 Normalized density profiles of oxygen atoms (red) and hydro-gen atoms (black) of water as functions of distance from β -AgI(0001) surfaces. Panels A and B are for ice (supercooledby 20 K) and liquid water (supercooled by 4 K), respectively,near the silver exposed face. Panel C is for liquid water (super-cooled by 20 K) near the iodide exposed face. On the distancescale, zero (vertical green line) coincides with the position ofthe centre of the outermost ion of the face. ρ0 is the averagedensity of oxygen or hydrogen atoms far from the surface.Note the shifted scale in Panels A and B. . . . . . . . . . . . . 50Figure 3.8 Normalized density profiles of oxygen atoms (red) and hy-drogen atoms (black) of liquid water as functions of distancefrom the Ag exposed β -AgI(0001) surface. Panel A is thezero charge case (water supercooled by 4 K). Panel B is forAgI carrying charges of ±1.0e (water supercooled by 20 K).The green vertical line at zero coincides with the position ofthe centre of the outermost ions, and ρ0 is the average densityof oxygen or hydrogen atoms far from the surface. Note theshifted scale for Panel A. . . . . . . . . . . . . . . . . . . . . 52xviiFigure 3.9 Snapshots of water on the γ-AgI(111) surfaces. Panel A is atop view of water (supercooled by 4 K) on the silver exposedface, and Panel C is a side view of a single hexagonal ring onthis surface. Panels B and D are corresponding views of water(supercooled by 20 K) on the iodide exposed face. Colours forsilver ions, iodide ions, oxygen atoms and hydrogen atomsare silver, green, red, and black, respectively. Yellow lineshighlight the hexagonal arrangement of the water on the surface 53Figure 3.10 Density profiles of oxygen (red) and hydrogen (black) atomsfor the six-site model near the γ-AgI(111) surface. Panels Aand B are for ice (supercooled by 20 K) and liquid water (su-percooled by 4 K), respectively, near the silver exposed face.Panel C is for liquid water (supercooled by 20 K) near the io-dide exposed face. On the distance scale, zero (vertical greenline) coincides with the position of the centre of the outermostion of the surface. ρo is the average density far from the sur-face. Note the shifted scale in Panels A and B. . . . . . . . . . 54xviiiFigure 3.11 Snapshots of liquid water on γ-AgI(001) surfaces. Panel Ais a top view of water (11 K above the melting point) on thesilver exposed face, and Panel C is a tilted view of a singleFCC-like arrangement on this surface. Panels B and D arecorresponding views of water (supercooled by 20 K) on theiodide exposed face. Panel E is a tilted view of a FCC ar-rangement in Ic. Silver ions, iodide ions, oxygen atoms, andhydrogen atoms are silver, green, red and black, respectively.The yellow lines highlight the FCC-like arrangement of thewater on the surface. . . . . . . . . . . . . . . . . . . . . . . 56Figure 3.12 Normalized density profiles of oxygen atoms (red) and hy-drogen atoms (black) of water as functions of distance fromγ-AgI(001)) surfaces. Panels A and B are for ice (supercooledby 20 K) and liquid water (6 K above the melting point), re-spectively, near the silver exposed face. Panel C is for liquidwater (supercooled by 20 K) near the iodide exposed face. Onthe distance scale, zero (vertical green line) coincides with theposition of the centre of the outermost ion of the face. ρ0 isthe average density of oxygen or hydrogen atoms far from thesurface. Note the shifted scale in Panels A and B. . . . . . . . 57xixFigure 3.13 Normalized density profiles of oxygen atoms (red) and hydro-gen atoms (black) of liquid water (supercooled by 20 K) asfunctions of distance from the Ag exposed γ-AgI(001) sur-face. Panel A is the zero charge case. Panel B is for AgIcarrying charges of ±0.8e. The green vertical line at zero co-incides with the position of the centre of the outermost ions,and ρ0 is the average density of oxygen or hydrogen atoms farfrom the surface. Note the shifted scale for Panel A. . . . . . . 58Figure 3.14 Snapshot of TIP4P/Ice water near the Ag exposed γ-AgI(001)face with charges of ± 1.0e on the AgI (supercooled by 20K). Colours for silver ions, iodide ions, oxygen atoms, andhydrogen atoms are silver, green, red, and black, respectively. . 59Figure 3.15 The first layer of water on the Ag+I exposed β -AgI(101¯0) face(supercooled by 20 K). Colours for silver ions, iodide ions,oxygen atoms, and hydrogen atoms are silver, green, red, andblack, respectively. The yellow and light blue lines highlightthe two orientations of the boat conformations made by thesilver iodide surface; yellow is opening outward and blue isopening inward towards the surface. . . . . . . . . . . . . . . 60xxFigure 3.16 Density profiles of oxygen (red) and hydrogen (black) atomsfor six-site water near the Ag+I exposed β -AgI(101¯0) surface(supercooled by 20 K). On the distance scale, zero (verticalgreen line) coincides with the position of the centre of theoutermost ion of the surface. ρo is the average density farfrom the surface. . . . . . . . . . . . . . . . . . . . . . . . . 61Figure 4.1 Top view of the AgI disks and plates considered. Silver andgreen spheres represent silver and iodide ions, respectively. . . 65Figure 4.2 Polarization profiles of water at ∆Ts values 1 K above whereice nucleation was observed. Red, green, and blue representthe X , Y and Z components, respectively. The system labeland the ∆Ts value are given in the legend. In all systems exceptC2, the silver exposed face of the AgI disk is directed outwardalong the positive Z axis. For C2, the silver exposed faces aredirected towards each other, such that any long-range fieldscancel. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69xxiFigure 4.3 Series of snapshots of an ice cluster nucleating bulk ice on AgIdisk D. At each time, side and top views of the ice cluster areshown. Silver and iodide ions are shown in silver and green,respectively. In the first and third columns, black, blue, andorange represent hydrogen atoms, the oxygen atoms of Ih, andthe oxygen atoms of Ic, respectively. In the second and fourthcolumns, yellow and purple represent the oxygen atoms of iceinside and outside the cylindrical volume defined by the diskradius, respectively. Liquid water molecules are not shown. . . 71xxiiFigure 4.4 ∆Ts required to nucleate ice in 20–200 ns versus the radiusof the nucleating disk or plate compared with previous simu-lation results and with the GT equation based on the proper-ties of real water. Heterogeneous nucleation results for AgIdisks or plates are shown in black (see legend) for the six-siteand TIP4P/2005 models. The solid black line indicates thelowest supercooling required to nucleate ice on infinite AgIslabs, and the dotted black curve represents a fit of our six-sitemodel results to the GT equation. Red downward trianglesand the red line show heterogeneous nucleation results[81]for the mW model on graphene disks and infinite slabs, re-spectively. Blue triangles are homogeneous nucleation resultsfor the TIP4P-E model obtained employing infinite cylindri-cal ice clusters.[66] The GT equation for cylinders[136] is alsoshown (dashed blue curve). Homogeneous nucleation resultsemploying spherical ice clusters are in green (see legend) forthe TIP4P/2005,[2] TIP4P/Ice,[2] and TIP4P[64] models. Thereal-water GT equation[136] for spherical nuclei is shown asa green dashed curve. . . . . . . . . . . . . . . . . . . . . . . 73Figure 4.5 Ice cluster size vs time for simulations done with neutral Agand I atoms (left panel) and with pressure coupling maintain-ing the pressure of the simulation at 1 bar. Keys show thesystem and ∆Ts for each simulation. . . . . . . . . . . . . . . 74xxiiiFigure 4.6 Caption on next page. . . . . . . . . . . . . . . . . . . . . . . 77Figure 4.6 Plots showing the number of water molecules identified as ice-like growing on each disk or plate at the minimum ∆Ts, whereice nucleation was observed, and at a temperature one degreewarmer. The ∆Ts values are given in the legends. Our es-timates of the minimum and maximum limits on the criticalcluster size (see text) are shown as horizontal lines and opencircles, respectively. . . . . . . . . . . . . . . . . . . . . . . . 78Figure 4.7 Ice cluster size versus time for the C2 case (six-site model) atthe ∆Ts values shown in the legend. . . . . . . . . . . . . . . 78Figure 4.8 Ice cluster size versus time for the TIP4P/2005 model. Theleft panel is for disk D and the right panel for disk E. The hor-izontal lines indicate our lower limit estimates of the criticalcluster size. Estimates of the upper size limit are shown ascircles. ∆Ts values are given in the legend. . . . . . . . . . . . 79Figure 4.9 Number of ice-like molecules as a function of time for fourdisks and plates. Curves are included for the total numberof ice molecules (black) and for the numbers inside (red) andoutside (green) the cylindrical volume defined by the disk orplate radius. In each panel, a circle denotes the point identifiedas the upper limit on the critical cluster size. . . . . . . . . . . 81xxivFigure 4.10 Ice cluster size versus time for the six-site model. The totalcluster size is represented by a black curve, the amounts ofice inside and outside a cylindrical volume defined by the diskradius are red and green, respectively. The system label andthe ∆Ts value are given in the legend. . . . . . . . . . . . . . . 82Figure 4.11 As in Figure 4.9, but for the TIP4P/2005 model. . . . . . . . . 83Figure 4.12 Caption on next page. . . . . . . . . . . . . . . . . . . . . . . 84xxvFigure 4.12 Comparison of critical cluster results from our simulationswith previous experimental and theoretical work. (A) Com-parison with previous simulations. Present results are shownas bars joining the upper and lower estimates for the criticalcluster and are labeled indicating the disk or plate employed.Results for the six-site and TIP4P/2005 models are shown inblack and pink, respectively. The green dotted curve repre-sents the real-water GT equation.[136] Green points (see leg-end) are homogeneous nucleation results for the TIP4P/2005,[2]TIP4P/Ice,[2] and TIP4P[64] models. Blue triangles are ho-mogeneous nucleation results for the six-site model obtainedin an applied electric field.[77] (B) Comparison with experi-mental estimates obtained using microemulsions[139] for ho-mogeneous nucleation (purple squares, see legend) and het-erogeneous nucleation (orange downward triangles). The openpurple squares are homogeneous nucleation results obtainedemploying levitated droplets.[140] Other symbols are as inpanel A. (C) Simulation (present work) and experimental re-sults plotted versus the scaled temperature τ . The symbolsand curves are as in panels A and B. . . . . . . . . . . . . . . 85Figure 4.13 An example of the peaks selected for the failed cluster plot(Figure 7). This simulation was for the six-site model withplate C at ∆Ts = 19 K. . . . . . . . . . . . . . . . . . . . . . 89xxviFigure 4.14 Sizes of subcritical clusters (six-site model) that failed to nu-cleate for different disks and plates as functions of ∆Ts. Crit-ical cluster size ranges for the smallest ∆Ts where nucleationoccurred are shown for each disk or plate as bars joining esti-mates of upper and lower limits. . . . . . . . . . . . . . . . . 90Figure 5.1 Number of ice molecules as a function of time for three differ-ent Al-surfaces, as indicated in the legends. Results are shownfor the TIP4P/Ice (top) and six-site (bottom) water models.The rigid-kaolinite, free-H, and free-OH surfaces are describedin the text. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96Figure 5.2 Configurational snapshot of ice grown on the Al-surface (free-H) with the six-site model. The simulation cell is orientedsuch that one is looking along the c-axis (into the basal plane)of ice Ih. Water molecules visible within the hexagons formedby the ice Ih lattice are from a few sheets of ice Ic. Mirroredkaolinite slabs are at the top and bottom of the simulation cell.Hydrogen, aluminum, silicon, and oxygen atoms are black,pink, yellow, and red, respectively, except for oxygen atomsreported as ice, which are blue. . . . . . . . . . . . . . . . . . 97xxviiFigure 5.3 Number of ice molecules as a function of time on Al-surfacesfor simulations begun under free-H conditions and then switchedto free-OH. The parent free-H parts are shown as gray curves,and the free-OH continuations are indicated in colour. Resultsare shown for the TIP4P/Ice (top) and six-site (bottom) watermodels. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99Figure 5.4 Caption on next page. . . . . . . . . . . . . . . . . . . . . . . 100Figure 5.4 Snapshots showing the beginning and end of each of the free-OH simulations begun with configurations taken from free-H simulations for the six-site model. Each pair of snapshotsshows the beginning (top image) and end (bottom image) ofthe free-OH simulation. The first time given above each pairof snapshots is the length of the initial free-H simulation, andthe second time is the length of the simulation under free-OHconditions. The colours for oxygen, hydrogen, aluminum, andsilicon atoms are red, black, purple, and yellow, respectively. . 101Figure 5.5 Snapshots showing the beginning and end of each of the free-OH simulations begun with configurations taken from free-Hsimulations for the TIP4P/Ice model. Format and colours areas in Figure 5.4. . . . . . . . . . . . . . . . . . . . . . . . . . 101xxviiiFigure 5.6 Four panels illustrating the behaviour of TIP4P/Ice water aswell as the kaolinite hydroxyl groups near Al-surfaces. Thesurfaces included are rigid-kaolinite (A), free-OH (B), free-Hbefore ice nucleation (C), and free-H after ice nucleation andgrowth (D). Each panel contains density profiles (top left), aconfigurational snapshot showing water and surface hydroxylgroups (top right), and a snapshot of the first layer of wa-ter near the Al-surface and the surface hydroxyl groups (bot-tom). Particular atoms are indicated with the same colourin all density profiles and snapshots: water hydrogen atoms(black), water oxygen atoms (red), hydroxyl hydrogen atoms(blue), and hydroxyl oxygen atoms (green). Density profilesare aligned with the snapshots to their right, and the zero onthe y-axis coincides with the crystallographic position of thehydroxyl oxygen atoms of the Al-surface. Density profilesare shown only for atoms that can move during the simula-tion. Therefore, panel A includes only density profiles for theoxygen and hydrogen atoms of water, and the density profileof the kaolinite hydroxyl oxygen atoms is shown only in panelB. The yellow hexagons in panels C and D (bottom) highlightexamples of hexagonal rings formed by water in the surfacelayer. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103xxixFigure 5.7 Four panels illustrating the behaviour of six-site water as wellas the kaolinite hydroxyl groups near Al-surfaces. The sur-faces included are rigid-kaolinite (A), free-OH (B), free-H be-fore ice nucleation (C), and free-H after ice nucleation andgrowth (D). Each panel contains density profiles (top left), aconfigurational snapshot showing water and surface hydroxylgroups (top right), and a snapshot of the first layer of wa-ter near the Al-surface and the surface hydroxyl groups (bot-tom). Particular atoms are indicated with the same colourin all density profiles and snapshots: water hydrogen atoms(black), water oxygen atoms (red), hydroxyl hydrogen atoms(blue), and hydroxyl oxygen atoms (green). Density profilesare aligned with the snapshots to their right, and the zero onthe y-axis coincides with the crystallographic position of thehydroxyl oxygen atoms of the Al-surface. Density profilesare shown only for atoms that can move during the simula-tion. Therefore, Panel A only includes density profiles for theoxygen and hydrogen atoms of water, and the density profileof the kaolinite hydroxyl oxygen atoms is shown only in PanelB. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104xxxFigure 5.8 Number of ice Ih and ice Ic molecules as a function of timegrowing from the Al-surface under free-H conditions. Resultsare shown for the six-site (240 K) and TIP4P/Ice (230 K) wa-ter models, as indicated in the legend. Note that ice Ih is byfar the dominant ice form. . . . . . . . . . . . . . . . . . . . 106Figure 5.9 Hydroxyl group configurations used in the prepared rigid sur-face simulations. Oxygen atoms are red, and hydrogen atomsare blue; only hydrogen atoms directed upward are shown. . . 109Figure 5.10 Additional ice- and liquid-prepared surfaces obtained for theTIP4P/Ice model. Hydroxyl oxygen atoms are red, hydroxylhydrogen atoms are blue, and only hydrogen atoms directedupward are shown. . . . . . . . . . . . . . . . . . . . . . . . 110Figure 5.11 Additional ice- and liquid-prepared surfaces obtained for thesix-site model. Hydroxyl oxygen atoms are red, hydroxyl hy-drogen atoms are blue, and only hydrogen atoms directed up-ward are shown. . . . . . . . . . . . . . . . . . . . . . . . . . 111Figure 5.12 Caption on next page. . . . . . . . . . . . . . . . . . . . . . . 114Figure 5.12 Si-surface without (A) and with (B) a layer of water. The oxy-gen atoms of water are blue, and hydrogen atoms are black.The oxygen atoms of kaolinite are red, and silicon atoms areyellow. The yellow and green lines in panel B outline theeight- and six-membered rings, respectively, formed by wateron the Si-surface. . . . . . . . . . . . . . . . . . . . . . . . . 115xxxiFigure 5.13 Configurational snapshots of ice nucleating and growing onthe Si-surface. The results shown are for the six-site modelat 240 K. The colours are as in Figure 5.2 except that oxygenatoms of ice Ih and ice Ic are blue and green, respectively. . . . 115Figure 5.14 Close-up of the snapshot at 700 ns from Figure 5.13. The redarrows and indices indicate planes in ice Ih, and the blue arrowtraces the layer of ice in contact with the Si-surface. . . . . . . 116Figure 5.15 A photograph of a model of the ice formed on the Si-surface.Black balls represent water oxygen atoms. The Si-surface isrepresented by the black background at the bottom of the picture.117Figure 5.16 Snapshots of ice nucleating (six-site model at 250 K) on thelarge Si-surface. Hydrogen, aluminum, silicon, and oxygenatoms are black, pink, yellow, and red, respectively, exceptfor oxygen atoms of ice, which are blue for ice Ih and greenfor ice Ic. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119Figure 6.1 Images of potassium feldspar surfaces considered. Silicon,aluminium, potassium and bridging oxygen atoms are yellow,teal, pink, and red respectively. Hydroxyl groups, and surfaceoxygen atoms in the +OH and +O surfaces are purple. The toprow displays the (001) surface and the bottom row displaysthe (010) surface. The type of surface is indicated at the topof each column. . . . . . . . . . . . . . . . . . . . . . . . . . 123xxxiiFigure 6.2 Density profiles of liquid water near feldspar at 240 K. Wateroxygen and hydrogen profiles are red and black, respectively.Hydrogen profiles of feldspar hydroxyl groups are green. Foreach panel, the upper profiles are for six-site water and thelower profiles are for TIP4P/Ice water. Note the density scalefor six-site water is shifted 200 particles/nm3 for clarity. Aboveeach panel the surface in question is indicated. . . . . . . . . . 126Figure 6.3 Plots of the number of ice molecules for several gibbsite sim-ulations. Panel A is for free-H simulations at 240 K withboth the six-site (dotted lines), and TIP4P/Ice (solid lines) re-sults included. Red indicates simulations on gibbsite, blue isfor ice-scaled gibbsite, and green is for kaolinite-scaled gibb-site. Most lines are near the bottom in panel A because noice nucleated, making it difficult distinguish them. Panel Bis for several ice-prepared simulations. Dotted lines indicatekaolinite-scaled gibbsite and solid lines indicate normal gibb-site. Blue lines are for TIP4P/Ice and red is for six-site water. . 128Figure 6.4 Ice growing on the kaolinite-scaled gibbsite(001) surface withthe TIP4P/Ice water model at 240 K. Oxygen, hydrogen, andaluminium atoms are red, black, and pink, respectively. Oxy-gen atoms of ice Ih determined by the CHILL algorithm areblue. The simulation cell is oriented to look down the c-axisof ice Ih. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129xxxiiiFigure 6.5 Density profiles of liquid water near gibbsite at 240 K. Wateroxygen, water hydrogen, and hydroxyl hydrogen profiles arein red, black, and green respectively. The upper profiles arefrom six-site water and the bottom profiles are from TIP4P/Icewater. Note the density scale for six-site water has been shifted200 particles/nm3 for clarity. From left to right the plots arefor gibbsite, ice-scaled gibbsite, and kaolinite-scaled gibbsite.Zero on the Z-axis coincides with the centre of the surface’shydroxyl oxygen atoms. . . . . . . . . . . . . . . . . . . . . 130Figure 6.6 Density profiles of TIP4P/Ice water (panel A), and six-site wa-ter (panel B), near ice-prepared surfaces at 240 K. Upper pro-files (note the density scale shifted by 200 particles/nm3 forclarity) are for gibbsite, and lower profiles are for kaolinite-scaled gibbsite. Red and black represent oxygen and hydro-gen density profiles, respectively. Zero on the Z-axis coin-cides with the centre of the surface’s hydroxyl oxygen atoms. . 132xxxivFigure 7.1 Weblogo[149] plot of the probabilities of residues in the P.syringae INP for residues 208-959. Position is the residuelocation in a repeated sixteen residue segment. Letters rep-resent amino acids, and the size of the letter shows the prob-ability at which that amino acid appears at that location inthe repeated segment. Blue indicates hydrophillic residues,green is neutral, and black is hydrophobic. Residues appear-ing here are: A=Alanine, D=Aspartic acid, E=Glutamic acid,F=Phenylalanine, G=Glycine, H=Histidine, I=Isoleucine, K=Lysine,L=Leucine, M=Methionine, N=Asparagine, Q=Glutamine, R=Arginine,S=Serine, T=Threonine, V=Valine, W=Tryptophan, Y=Tyrosine.136Figure 7.2 Sequences and images of the proteins used in this study. PanelA displays the sequence (top panel) of the original AFP pro-tein from L. perenne, and an image (bottom panel) of the pro-tein. Residues in the sequence shown in red were removedduring mutation, and the residues in blue were changed. PanelB displays the sequence and image of the mutated protein usedin the simulations. Residues shown in blue were changed fromthe original AFP. In the image in Panel B threonine residuesare red, valine residues are blue, and serine residues are orange. 138xxxvFigure 7.3 Two views of the raft of proteins. Panel A shows a top viewand panel B shows an end-on view. The proteins are situatedin the simulation cell to produce an infinitely periodic surfaceof proteins. . . . . . . . . . . . . . . . . . . . . . . . . . . . 140Figure 7.4 Two orientations of the protein used in an attempt to “force”the protein into an icelike configuration. The protein (teal) isoriented either approximately parallel to the AgI disk (greenand silver) as in panel A, or perpendicular to the disk as inpanel B. Water molecules are not shown. . . . . . . . . . . . . 141Figure 7.5 Horizontal (Panel A) and vertical (Panel B) “ice-prepared”proteins. Colours of the proteins are the same as in Figure7.2, and the structure of threonine, valine, and serine residuesis displayed. Pink shapes show the regions occupied by wateroxygen atoms. . . . . . . . . . . . . . . . . . . . . . . . . . . 142Figure A.1 The number of ice molecules as a function of time for theTIP4P/Ice water model detected by the CHILL algorithm insimulations on three AgI surfaces: Ag exposed β -AgI(0001)(red), Ag exposed γ-AgI(001)( blue), and Ag exposed γ-AgI(111)(green). The solid lines and dotted lines represent the numberof Ic and Ih molecules, respectively. . . . . . . . . . . . . . . 165xxxviFigure A.2 Density profiles of oxygen (red) and hydrogen (black) atomsfor the TIP4P/Ice model near β -AgI(0001) surfaces. PanelsA and B are for ice (supercooled by 20 K) and liquid water(5 K above the melting point), respectively, near the silverexposed face. Panel C is for liquid water (supercooled by 20K) near the iodide exposed face. On the distance scale, zero(vertical green line) coincides with the position of the centreof the outermost ion of the surface. ρo is the average densityfar from the surface. Note the shifted scale in Panels A and B. 166Figure A.3 Density profiles of oxygen (red) and hydrogen (black) atomsfor the TIP4P/Ice model near γ-AgI(001) surfaces. Panels Aand B are for ice (supercooled by 20 K) and liquid water (10 Kabove the melting point), respectively, near the silver exposedface. Panel C is for liquid water (supercooled by 20 K) nearthe iodide exposed face. On the distance scale, zero (verticalgreen line) coincides with the position of the centre of theoutermost ion of the surface. ρo is the average density farfrom the surface. Note the shifted scale in Panels A and B. . . 167xxxviiFigure A.4 Density profiles of oxygen (red) and hydrogen (black) atomsfor the TIP4P/Ice model near γ-AgI(111) surfaces. Panels Aand B are for ice (supercooled by 20 K) and liquid water (5 Kabove the melting point), respectively, near the silver exposedface. Panel C is for liquid water (supercooled by 20 K) nearthe iodide exposed face. On the distance scale, zero (verticalgreen line) coincides with the position of the centre of theoutermost ion of the surface. ρo is the average density farfrom the surface. Note the shifted scale in Panels A and B. . . 168Figure A.5 Density profiles of oxygen (red) and hydrogen (black) atomsfor the TIP4P/Ice model near the Ag+I exposed β -AgI(101¯0)surface (supercooled by 20 K). On the distance scale, zero(vertical green line) coincides with the position of the centreof the outermost ion of the surface. ρo is the average densityfar from the surface. . . . . . . . . . . . . . . . . . . . . . . . 169Figure A.6 Density profiles of oxygen (red) and hydrogen (black) atomsfor the TIP4P/Ice model near the Ag exposed β -AgI(0001)surface . Panel A is for neutral AgI (5 K above the meltingpoint) and Panel B is for AgI with charges of ± 1.0e (super-cooled by 20 K). On the distance scale, zero (vertical greenline) coincides with the position of the centre of the outermostion of the surface. ρo is the average density far from the surface.170xxxviiiFigure A.7 Density profiles of oxygen (red) and hydrogen (black) atomsfor the TIP4P/Ice model near the Ag exposed γ-AgI(001) sur-face (supercooled by 20 K). Panel A is for neutral AgI andPanel B is for AgI with charges of ± 1.0e. On the distancescale, zero (vertical green line) coincides with the position ofthe centre of the outermost ion of the surface. ρo is the averagedensity far from the surface. . . . . . . . . . . . . . . . . . . 171Figure B.1 Orientations of water models tested to produce interaction en-ergies. Purple circles represent the ion, either Ag or I, and redand black circles represent the oxygen and hydrogen atoms,respectively, of water. Orientation A has the dipole momentof water pointing toward the ion. In B, the dipole momentof water is pointing away from the ion, and in C a hydrogenbond is pointing toward the ion. The arrow labeled ‘r’ showsthe direction the water was moved in the calculation of theinteraction energies. . . . . . . . . . . . . . . . . . . . . . . . 179xxxixFigure B.2 Interaction energies between an Ag or I ion and a six-site wa-ter model in the three orientations shown is Figure B.1. Thetitle of each plot lists the orientation and the ion involved. ZBPcurves are from the potential used in this study. ZBP 0.8e andZBP 0.0e are the present potential using charges of ± 0.8eand 0.0 on the AgI ions, respectively. HK is the potential[35]from which the ZBP potential is derived from and includespolarization terms. . . . . . . . . . . . . . . . . . . . . . . . 180xlList of Symbols and Abbreviationsa Accelerationai Lattice parameter of icean Lattice parameter of an ice nucleusaq CHILL order parameterb Simulation cell geometryb1/b2/b3 Simulation cell edge vectors∆c Difference between isobaric specific heats of water and iceC Nucleation rate pre-exponential factorC1/C2 Constants for Gibbs-Thompson equationd Degree of B-spline interpolationdOH Distance between oxygen and hydrogen atomsdOL Distance between oxygen atom and lone pairdOm Distance between oxygen atom and m-sitee Elementary chargef Degrees of freedomxlif (mθ ) Heterogeneous nucleation form factorF(Q) Fourier transform of charge arrayF Forceg Gibbs–Thomson geometry parameter∆g∗ Activation energy of monomer attaching to a cluster∆G Change in free energyIc Cubic iceIh Hexagonal iceIsd Stacking disordered iceIm Lattice mismatch integerIn Lattice mismatch integerJ Nucleation ratek SPME mesh pointK Number of SPME mesh pointskB Boltzmann constantL Lagrangianmi Particle massm(z) Average dipole order parametermθ Cosine of contact angleMd(η) B-splinen Simulation periodic imagesn1/n2/n3 Integers representing periodic imagesni Number density of icexliiN Number of particles in a simulationNl Density of waterp MomentumP PressurePre f Reference pressureq Heat of meltingqH Hydrogen atom partial chargeqi/q j Partial chargeqlm CHILL local order parameterqLP Lone pair partial chargeqm m–site partial chargeQ Charge arrayQm Nose´-Hoover mass parameterr Ice embryo radiusr i Particle coordinatesr i j Inter-particle distances Nose´-Hoover virtual particle coordinates j Fractional coordinateS(w) Structure factorS˜(w) Approximate structure factort Time∆t Time between time stepsT TemperaturexliiiTm Melting temperature∆Ts SupercoolingTw Temperature of the Widom lineui j Potential energyuCoul Coulombic potential energyuIND Polarization potential energyuLJ Lennard-Jones potential energyulong Long-range Coulombic potential energyurec Reciprocal Coulombic potential energyusel f Coulombic self potential energyushort Short-range Coulombic potential energyvi VelocityV Simulation cell volumew Reciprocal lattice vectorsw1/w2/w3 Integers for reciprocal lattice vector imagesW Strength of pressure couplingYl,m Spherical harmonicα Gaussian distribution width factorαm Polarizability of AgI ionsαw Polarizability of waterγi,n Surface tension between ice and an ice nucleusγl,i Surface tension between ice and liquid waterγl,n Surface tension between liquid water and an ice nucleusxlivδ Lattice mismatchε0 Vacuum permittivityεH Hydrogen atom energy well depthεi j Lennard-Jones energy well depthεO Oxygen atom LJ energy well depthεi/ε j Lennard-Jones particle energy well depthζ Thermodynamic friction coefficientη Fractional scaled coordinateθ Contact angleκ Thermal expansivity of iceµ i Dipole moment vector∆µl,i Change in chemical potential between liquid and iceµw Dipole moment of waterρ Densityρ0 Bulk water densityσH Hydrogen atom radiusσi j Lennard-Jones radiusσi/σ j Lennard-Jones particle radiusσO Oxygen atom radiusτ Scaled temperatureχ Compressibility of iceAFP Anti-freeze proteinCCN Cloud condensation nucleixlvCNT Classical nucleation theoryCPU Central processing unitFCC Face centred cubicGT Gibbs-Thompson equationHK Hale and Kiefer AgI force fieldIN Ice nucleusINP Ice nucleating proteinLJ Lennard-JonesMD Molecular dynamicsNPT Constant number of particles, pressure, and temperatureNVT Constant number of particles, volume, and temperatureSPME Smooth particle mesh EwaldZBP Zielke-Bertram-Patey AgI force fieldxlvi4 years real time∼463 years, 337 days, 4 hours, 22 minutes processor time∼148.775 µs simulated< /PhD >xlviiChapter 1Introduction“Now suppose,” chortled Dr. Breed, enjoying himself, “that therewere many possible ways in which water could crystallize, couldfreeze. Suppose that the sort of ice we skate upon and put intohighballs–what we might call ice-one–is only one of several types ofice. Suppose water always froze as ice-one on Earth because it hadnever had a seed to teach it how to form ice-two, ice-three, ice-four...?And suppose,” he rapped on his desk with his old hand again, “thatthere were one form, which we will call ice-nine–a crystal as hard asthis desk–with a melting point of, let us say, one-hundred degreesFahrenheit, or, better still, a melting point of one-hundred-and-thirtydegrees.” — Kurt Vonnegut, Cat’s CradleWhy is the sky blue? What came first, the chicken or the egg? Do you knowabout Tyler Durden?To this list of classic questions I would like to add: How does water freeze?I add this question because a cursory examination of the relevant kinetics andthermodynamics of nucleation reveals that the freezing of water is a much moredifficult process than is commonly believed. Left to its own devices, water will not1freeze homogeneously until very cold temperatures. For example, small dropletswill not freeze homogeneously until−38◦C.[1] As another example, one estimatefrom simulations suggests that if all the water on earth was given the age of theuniverse to freeze homogeneously, then -20◦C is the warmest one could expect iceto nucleate.[2]1.1 Kinetics of Ice NucleationClassical nucleation theory (CNT) is a framework commonly used to describenucleation. In CNT for homogeneous nucleation, it is usually assumed the nu-cleating phase is spherical. In the case of water, it is assumed, a spherical iceembryo nucleates out of the surrounding liquid water, producing a change in thefree energy of the system given by eq 1.1[3, 4]∆G(r)homo = ni4pir33(∆µl,i)+4pir2γl,i, (1.1)where the change in free energy of the system, ∆G, is a function of the radius,r, of a spherical ice embryo forming in liquid water. In eq 1.1, ni is the numberdensity of ice, ∆µl,i is the change in chemical potential for the bulk liquid to icephase transition at a particular temperature, and γl,i is the surface tension betweenice and water. The first term on the right-hand-side of eq 1.1 produces a negativevalue, driving the spontaneity of the process. The second term produces a positivevalue, destabilizing the growth of the embryo.Figure 1.1 is a plot of eq 1.1 using data from Espinosa et al.[5] gathered from2molecular dynamics simulations,[6] demonstrating the behaviour of the equation.The figure shows that at small radii, the second term in eq 1.1 (surface tensionterm) is larger (because r2 > r3 for small numbers). However, as the radius be-comes larger the first term dominates. Figure 1.1 (black line) shows that initiallythe change in free energy increases as the radius of the ice embryo increases untilthe curve reaches a maximum, referred to as ∆G∗homo. The embryo at the maxi-mum is referred to as the critical cluster. Additional growth beyond the criticalcluster decreases the change in free energy, hence the critical cluster will tend togrow into bulk ice. In practice, cold temperatures are required to nucleate ice be-cause the change in chemical potential increases with decreasing temperature, anda large change in chemical potential is required to overcome the surface tensionof small embryos.Figure 1.1: Plot of eq 1.1 (black line). The red line is the surface energyterm and the blue line is the absolute value of the chemical potentialterm.The rate of critical ice clusters forming homogeneously per volume of liquid3water is given by,J =Cexp−∆g∗kBT Nlexp−∆G∗homokBT , (1.2)where J is the nucleation rate in units of m−3s−1.[7] C is a constant defined by therate of consumption of clusters as embryos grow. Nl is the density of liquid wa-ter and the second exponential factor accounts for the thermodynamic parameterspreviously introduced. The first exponent is the classical transition frequency[4]and describes the kinetic portion of the rate from the rate of attachment of watermolecules to the embryo. The form of the first exponent is a Boltzmann distri-bution, where ∆g∗ is the activation energy of a water molecule attaching to thecluster.[7] ∆g∗ is a function of the number of hydrogen bonds a water moleculeis participating in, and therefore varies with temperature, getting larger at lowertemperatures. Eq 1.2 gives the rate of nucleation as a product of the number ofcritical clusters in a volume of water, Nl multiplied by the second exponential, andthe rate of attachment of water molecules to an existing cluster, first exponential.Water can freeze at temperatures warmer than −38◦C with the aid of a het-erogeneous ice nucleus (IN). The IN acts as a catalyst, lowering the change infree energy of ice embryo formation by providing a more favourable interfacewith the ice embryo than liquid water. Mathematically, the effect of the IN canbe represented by a factor, f (mθ ), multiplied by the change in free energy for thehomogeneous formation of an embryo to produce ∆Ghet , the change in free energyfor heterogeneous formation of an embryo,4∆Ghet(r) = ∆Ghomo(r) f (mθ ), (1.3)mθ = cosθ =γln− γinγli, (1.4)where mθ ranges from -1 to 1, and is a function of the contact angle, θ . Thecontact angle is the angle produced between an ice embryo and an IN surface. mθcan also be calculated from the surface tensions between the three componentsof the system: liquid water and ice, γli; liquid water and the IN, γln; ice and theIN, γin.[8, 9] The form of the function f (mθ ) depends on the geometry of thenucleating system, e.g, f (mθ ) = (2 + mθ )(1−mθ )2/4 for a hemisphere of icenucleating on a flat surface.[10]1.2 Atmospheric Aerosols as Ice NucleiAirborne particles, or atmospheric aerosols, are a vital yet difficult to understandcomponent of Earth’s climate and weather. Atmospheric aerosols include a largenumber of different particle types, including: bacteria, fungal spores, mineral dustparticles, soot, sodium chloride/sea salt, sulfates, and nitrates. Aerosol sizes rangefrom 1 nm to 10 µm or more. The number and type of aerosols in a given areadepend on the emission sources, time of year or day, and weather conditions.[1, 3]Climate can be directly affected by aerosols which can scatter or absorb sun-light. Absorbed light will lead to heating and scattered light results in cooling,since a portion of scattered light is directed back into space.[3] The net effect of5aerosols on climate due to absorbing and scattering light is still uncertain.[11]Aerosols also affect climate indirectly by nucleating clouds. Clouds can reflectsunlight away from earth, which has a cooling effect on climate, or trap terrestrialradiation, leading to heating. Pollution affects clouds by increasing the numberof cloud condensation nuclei (CCN) in the atmosphere. Higher numbers of CCNlead to clouds with larger numbers of smaller droplets, which are more reflective.Secondly, these smaller droplets are less likely to precipitate, thus the clouds lastlonger.Without the presence of CCN, clouds would not form until exceedingly highpartial pressures of gas-phase water.[3] CCN help cloud droplets to nucleate throughRaoult’s law, which states that a dissolved species in a liquid will decrease thevapour pressure of the liquid compared to the pure liquid. By dissolving in water,CCN lower the vapour pressure of the droplet. CCN are therefore often solublesalts or molecules, in contrast to IN, which are usually insoluble.Often formed with the aid of IN, ice plays an important role by producing amajority of rainfall via the Bergeron-Findeisen process.[3, 12] When water vapourin the air is subsaturated with respect to liquid water, but saturated with respectto ice, liquid water droplets will evaporate and ice particles will grow. Once largeenough, the ice particles will start to fall as precipitation.Particles that are known to be good heterogeneous ice nuclei include mineraldusts, certain biological particles, and silver iodide.[1, 13] However, why someparticles function as good ice nuclei while others do not is poorly understood.The absence of a good molecular level understanding of heterogeneous ice nucle-6ation means that researchers cannot predict the ice nucleation ability of a substratewithout experiments. This absence of a microscopic understanding can lead to un-certainties when modelling ice nucleation in the atmosphere.[14]Macroscopically, there are three mechanisms by which an aerosol can nucleateice:[3] deposition nucleation, water vapour freezes immediately upon condensingonto an IN; immersion freezing, an IN inside a liquid drop initiates freezing; con-tact freezing, an IN collides with a droplet, causing the droplet to freeze.A number of requirements have been suggested for an aerosol to be an effec-tive IN based on experimental observations.[1] First, the IN must be insoluble orelse it will dissolve and not be able to provide a surface for nucleation. Second, theIN size is thought to be important. Larger IN are usually more active than smallerIN and positive correlations have been observed between size and ice nucleatingefficiency. Third, interactions between the IN and water need to be similar instrength and behaviour to inter-water hydrogen bonds. Fourth, an IN should havea crystal lattice which matches that of ice. This match should be both in size andshape of the unit cell. Differences in lattice parameters can be calculated usingthe following equation,[1]δ =Inan− ImaiImai×100%, (1.5)where δ , is the mismatch between the ice unit cell, ai, and the unit cell of theIN, an. Im and In are integers chosen to minimize δ . Experimental evidence fromthe nucleation of a variety of solids suggests the smaller the mismatch between7a nucleation catalyst, the IN in the case of ice, and the solid phase the better thecatalyst nucleates the solid phase.[1, 15]Fifth, a requirement of a good IN is the presence of active sites. Active sitesare locations on an aerosol’s surface possessing some property which makes thesite especially effective at nucleating ice. One study on dust particles found threedifferent types of active sites; the most efficient active sites being the rarest.[16]What active sites are is currently unknown or if other IN requirements, such asmismatch, are important for an active site’s operation.1.3 IceDepending on temperature, pressure, and the conditions ice forms, ice can existin many different forms. The phase diagram of water currently contains twelvecrystalline phases of ice,[17] and three forms of amorphous ice.[17, 18] This countexcludes clathrates, of which there are several. In addition, recently a phase of icewhich appears only during nucleation has been proposed, and is called ice 0.[19]At tropospherically relevant temperatures and pressures (1-0 atm, 0◦C to -50◦C),ice I is the thermodynamically stable phase.Ice I can exist in two forms, hexagonal ice, Ih, and cubic ice, Ic. Ice Ih isthermodynamically more stable than the metastable ice Ic.[17, 18] If formed inexperiments, ice Ic converts to ice Ih over the course of minutes to hours.[20, 21]Both phases of ice I can be constructed by stacking bilayers of water formed ofhexagonal rings in the chair conformation. The phases differ in how the bilay-ers are stacked.[22] Ice Ih has ABAB... stacking resulting in a hexagonal unit8cell, whereas, ice Ic has ABCABC... stacking giving a FCC (face centred cubic)lattice.[18, 22] The crystal structures are such that the Ih(0001) and Ic(111) facesare almost identical and can connect with little lattice strain.[23] Because of thecommon bilayer, the occurrence of random layers of Ih and Ic in a single ice crys-tal is observed in both simulations[22] and experiments,[24, 25] and is referred toas stacking disordered ice, Isd .[26]Both ice Ih and Ic crystal structures are based only on the position of water oxy-gen atoms. Hydrogen atoms in ice are usually randomly arranged, which wouldleave some residual entropy even at 0 K. However, at very cold temperatures, IceXI can form, which is a proton ordered phase of ice Ih.[6]1.4 Ice Nucleating MaterialsSilver iodide is one of the best heterogeneous ice nucleating substances known,inducing ice formation at temperatures as high as−3◦C.[1] Silver iodide is some-times used as a cloud seeding agent.[27, 28] Early researchers often assumed thatthe good ice nucleation activity of silver iodide was due to the small lattice pa-rameter mismatch between hexagonal ice and silver iodide surfaces.[1, 29] Sub-sequent work, however, has called into question the assumption that a small latticemismatch is indicative of ice nucleation ability.[30, 31]Two polymorphs of silver iodide can exist under atmospheric conditions. Theβ phase is the more stable and has a hexagonal unit cell.[32, 33] The γ phaseis metastable and has a cubic unit cell.[32] Most theoretical studies of the icenucleation ability of silver iodide have focused on the β -AgI(0001) face due to its9hexagonal symmetry, and this face has been shown to order water into hexagonalrings.[34–43]Every year approximately 1000-4000 Tg of mineral dust is carried up into ouratmosphere.[11] This mineral dust is believed to be a large source of atmosphericheterogeneous ice nuclei,[13, 44] yet a good understanding of how dust particlesnucleate ice is lacking.[45] Experimental studies have revealed certain types ofmineral dust such as illite, feldspar, and kaolinite, to be effective ice nuclei.[13,44, 46] However, the microscopic mechanisms responsible for ice nucleation onthese surfaces have not been determined. Airborne particles of these dusts nodoubt contain a large variety of surface structures and active sites, with varyingice-nucleating ability.[16]Kaolinite, Al2Si2O5(OH)4, is a type of mineral dust and also a clay. Thismaterial has been reported to nucleate ice at temperatures up to −10◦C,[47, 48]although typical nucleation temperatures are below −20◦C.[13, 44] The primarycleavage plane of kaolinite is the (001) plane,[49] which has two different sur-faces. One surface, referred to here as the Al-surface, is covered by hydroxylgroups attached to Al atoms arranged in a hexagonal pattern. The other surface,referred to as the Si-surface, consists of silicon atoms and bridging oxygen atomsforming hexagonal rings. In simulations the Al-surface has been shown to behydrophilic and the Si-surface hydrophobic.[50, 51]Potassium Feldspar, KAlSi3O8, (also a type of mineral dust) has recently beenfound to be an extremely effective IN in experimental studies.[46, 52] Later stud-ies have shown that the ice nucleating efficiency of a dust sample increases with10increasing potassium feldspar content.[53]Gibbsite, Al(OH)3, is a weathering product of kaolinite, feldspar and otheraluminosilicates. Under certain conditions, precursor species loose their siliconatoms, leaving only aluminium and hydroxide ions.[54–58] The (001) face ofgibbsite bears remarkable similarities to the Al-surface of kaolinite suggestingit should have ice nucleating ability. We are not aware of any experimental resultson the ice nucleation ability of gibbsite.Some of the most effective natural IN are ice nucleating proteins created bycertain types of bacteria. The most commonly studied ice nucleating proteincomes from the bacteria Pseudomonas syringae.[59] No crystal structure of thisprotein has been reported, however the sequence of amino acids is known.[60]The sequence shows that most of the protein is composed of a repeating sequenceof amino acids, likely in some sort of a helical structure. The P. syringae proteinalso contains end groups, which do not share in the repeating structure of the in-terior of the protein, and are believed to help mount the protein on the bacteria’ssurface. Another study showed that clusters of these proteins were necessary tonucleate ice, with 50-60 capable of nucleating ice at -3◦C.[61] A computationalstudy using an estimated structure of the ice nucleation protein from PseudomonasBorealis,[62] whose sequence is very similar to ice nucleation protein from P. sy-ringae, showed that the sequence is very stable as a β -helix, a helix forming two,parallel, β -sheets.111.5 Previous Simulation Studies of Ice NucleationComputer simulations are rapidly emerging as a valuable tool in the ongoing effortto understand the microscopic mechanisms of ice nucleation, and nucleation pro-cesses in general.[63] Homogeneous nucleation has been investigated with atom-istic water models,[2, 64–67] and with the coarse grained (one site per molecule)mW water model.[6, 68] It is very difficult to nucleate ice homogeneously, espe-cially using atomistic water models. The mW model[6] is advantageous in homo-geneous nucleation studies because its dynamics are very fast. The mW model,based on the Stillinger-Weber potential for silicon,[69] contains only one site anduses a three body term to govern the position of neighbouring water molecules,forcing the molecules to adopt icelike positions. Being a single particle withoutcharges, means that the mW model is very computationally cheap to run.[6]Observing homogeneous ice nucleation in direct simulations with atomisticmodels is very difficult because the probability of nucleation is extremely lowat higher temperatures (but still below melting), and at colder temperatures thedynamics of water models decreases to such a point that very long simulationsare required to explore sufficient phase space to achieve nucleation. To get aroundthese problems with atomistic models, researchers have used methods such as cre-ating prefabricated ice clusters,[2, 64–66], forward flux sampling,[70], and bias-ing potentials.[71] It has been suggested that a simulation of 20 000 TIP4P/2005water models[72], a type of atomistic model, running for 6 µs at 195 K wouldfreeze homogeneously.[64] The scale and time of this simulation however is be-12yond the simulations in this thesis.Electric fields have been shown to nucleate ice in simulations,[73–77] how-ever, experimental evidence for field induced nucleation is conflicting,[78, 79]possibly due to the fact that the electric fields used in simulations are too largefor experiments thus smaller fields are used. Ice nucleated by an electric field insimulations is Ic because electric fields polarize ice[74] and Ic can exist in polarform, whereas Ih does not form under polarizing conditions.[76] However, if onlya portion of a simulation cell is exposed to the field, Ih may appear in ice growingin regions of the cell untouched by the field.[76] Applying an electric field has theeffect of raising the melting point of water,[80] allowing for the water to be deeply“supercooled” while still at a temperature warm enough to allow the dynamics tobe fast enough to observe homogeneous nucleation if the field covers the entiresample.[77]If a good IN can be found, then ice nucleation poses an easier task for di-rect simulations, and spontaneous nucleation followed by ice growth has beenobserved on infinite sheets and disks of graphene[81, 82] with the mW model,[6]and on infinite sheets of silver iodide with atomistic models (Chapter 3).[83] Inthe case of graphene, the authors explicitly remark[81] that the effect is not dueto templating of the ice structure by the carbon surface, but due to the formationof layers of water next to the surface. Cox et al.[84, 85] have reported simula-tions of ice nucleation on hemispherical nanoparticles and disks employing themW water model.[6] They found that the ability of a surface to organize watermolecules into an icelike configuration was controlled by the surface morphology13and the water-surface interaction energy. Optimization of the surface morphol-ogy and interaction energy can increase the nucleation rate significantly over thehomogeneous ice nucleation rate.Ice nucleation with atomistic models influenced by a kaolinite surface has alsobeen reported,[86] but attributed to the influence of order imposed by the surfacein the adjacent water layers, evident in water density profiles, rather than to anyparticular match between ice and the underlying surface. Additionally, the simu-lation results described in reference [86] strongly depended on system size, andice nucleation could not be achieved for larger systems. Previous Monte Carlosimulations suggest that a rigid Al-surface of kaolinite (001) is not capable of or-ganizing water into icelike configurations.[87, 88] Instead, it was suggested thattrenches or surface defects may be responsible for ice nucleation as they can pro-duce electric fields that order the water within them.[89, 90] Quantum mechanicsbased calculations have shown the Al-surface of kaolinite (001) to be amphotericas the hydroxyl groups rotate allowing water to both donate and receive hydrogenbonds from the surface.[51, 91–96] No bulk ice nucleation was observed in thesestudies, partially due to the small size of these simulations, although the formationof an icelike monolayer on the Al-surface was observed.[94]1.6 Objectives and OutlineExperimentally, a number of effective IN have been discovered along with thetemperatures at which the IN are active. What is missing from our understandingis a microscopic picture of ice nucleation on an IN. We would like to determine14why some IN are more effective than others.We use molecular dynamics simulations to probe heterogeneous ice nucle-ation by several atmospherically relevant materials and to improve our micro-scopic understanding of heterogeneous ice nucleation. General methods are givenin Chapter 2, while methods specific to each chapter are presented in the relevantchapters. First we explore silver iodide. Being one of the best IN known, an un-derstanding of how silver iodide nucleates ice may give us a basis to understandhow other materials nucleate ice. In Chapter 3 we explore seven infinite surfacesof silver iodide. Three surfaces nucleated ice by organizing water into an icelikeconfiguration. We then study finite disks and plates of silver iodide in Chapter 4and find that larger disks are capable of nucleating ice at warmer temperatures bysupporting larger critical clusters.Kaolinite is studied in Chapter 5. This work shows that in some cases boththe IN surface and ice are able to adapt themselves to facilitate ice nucleation.We then study feldspar and gibbsite in Chapter 6, exploring the possibility thatweathering products from feldspar may be a source of IN. Lastly, in Chapter 7 wereport preliminary studies on the ice nucleating ability of proteins.15Chapter 2General Methods“All of the true things I am about to tell you are shamelesslies.”— Kurt Vonnegut, Cat’s CradleWhile computers have gotten significantly faster in recent years, calculationsinvolving large numbers of atoms are still computationally expensive. Quan-tum mechanics and Density Functional Theory methods require large amountsof computer resources, limiting the size and time scales of simulations. Classi-cal molecular dynamics (MD), where interparticle interactions are composed ofsimple equations, allow for simulation of 105 or more particles for time scales ofhundreds to thousands of nanoseconds. Heterogeneous ice nucleation, as investi-gated in this thesis, occurs on time scales of tens to hundreds of nanoseconds, insystems containing thousands of molecules. For these reasons we use the GRO-MACS 4.5.5 and 4.6.2 software packages.[97] GROMACS is intended for MDsimulation of biological molecules (proteins, DNA/RNA, lipids) and is designedto handle simulations of the size we require. Using available software also saves16time in programming and testing a new code from scratch, especially for the highdegree of paralellization required (often 32 or 64 CPU’s per simulation).2.1 Molecular DynamicsTo create a molecular dynamics simulation using GROMACS, several files areread into a pre-processor program including a coordinate file, containing the coor-dinates of all the atoms in the simulation; topology file, which specifies the type,geometry, number, and interactions of the atoms and molecules; parameter file,which specifies the conditions of the simulation, length, temperature, constraints,etc...; index file, containing the positions in the coordinate file of various groupsof atoms. The pre-processor produces a run input file which is fed into the mainmolecular dynamics program.A MD simulation begins with initialization of the program.[98] Initializationloads the relevant data into the program. Initial velocities may be set either froma Maxwell-Boltzmann distribution or read in. For parallel simulations, certainCPUs may be dedicated to particular tasks.The main loop of a MD program calculates the forces between particles, thenmoves the particles according to the equations of motion. This step of the programrequires almost all of the time required to run a simulation, and significant workhas been done to make these calculations as efficient, yet accurate, as possible.[98]Calculation of forces and moving the particles is repeated until the end ofthe simulation. A trajectory file is produced during the simulation containing thepositions of the particles at regular intervals, and sometimes velocities and forces.172.2 InteractionsIn our simulations particles are either single site atoms or ions, or molecules madeup of several sites. Sites in a molecule can represent the atoms in the molecule,or artificial sites added to give the molecules more realistic behaviour. The in-teraction energy between two particles is a function of the distance between thesites within each particle. For the models we employ, the energy between sites isa combination of Coulombic and Lennard-Jones (LJ) interactionsui j(r i j) =qiq j4piε0r i j+4εi j[(σi jr i j)12−(σi jr i j)6], (2.1)where the interaction energy, ui j, between sites i and j is a function of the dis-tance between them, r i j.[98] The first term on the right-hand-side is the Coulom-bic energy from partial charges, q, on each site. ε0 is the vacuum permittivity. Thesecond term is the LJ potential which is repulsive at short distances to give the par-ticle size, and attractive over longer distances to mimic London dispersion forces.σi j determines the “diameter” of the site, and εi j sets the depth of the potentialenergy well. A force field determines the charges and values of σi and εi for eachsite. Not every site in a molecule contains values for both Coulombic and LJ inter-actions. Some sites are point charges, which only have the Coulombic term, andothers only have LJ parameters. When two particles with different values of σiand εi interact, the values used in eq 2.1 are calculated using the following mixingrules,[99]18σi j =σi +σ j2, and εi j =√εiε j. (2.2)σi j is an arithmetic mean of the values from two interacting sites, and εi j is ageometric mean.When two particles are interacting, the total energy between them is the sum ofthe interactions between all of the sites in one particle and all the sites in the otherparticle. As an example, consider the TIP4P/Ice[100] water model used in thisthesis (Figure 2.1). Additional details of this model are given later in this chapter.Briefly, the model is comprised of four sites: an oxygen atom with only a LJ term,two hydrogen atoms which are positive point charges, and an m-site which is anegative point charge. The energy between two TIP4P/Ice water models A and Bisu(A,B) = uLJ(AO,BO)+uCoul(AH1,BH1)+uCoul(AH1,BH2)+uCoul(AH1,Bm)+uCoul(AH2,BH1)+uCoul(AH2,BH2)+uCoul(AH2,Bm)+uCoul(Am,BH1)+uCoul(Am,BH2)+uCoul(Am,Bm).(2.3)Simulations are conducted in a simulation cell whose geometry is denoted bythe matrix b which has vectors b1,b2,and b3, which define the edges of the cell.This cell is repeated in every direction an infinite number of times producing aninfinite periodic system. To avoid having to calculate energies between all parti-19Six-site TIP4P/IceFigure 2.1: Six-site and TIP4P/Ice water models. Oxygen atoms, hydrogenatoms, lone pairs, and m-site are red, black, teal, and pink respectively.TIP4P/2005 shares the same basic shape as TIP4P/Ice.cles and their images, the LJ term in eq 2.1 is truncated at 1 nm. The LJ potentialdecays rapidly, 1/r6, which is faster than the volume increases with increasingdistance,[98] ∝ r3, so the contribution of the LJ potential past the cut-off is negli-gible. The value of the LJ potential at the cut-off is added to the potential over thepreceding distance so the LJ potential is zero at the cut-off; this is known as thecut and shift method.With periodic boundary conditions, the Coulombic contribution to the poten-tial energy becomes[98, 101]uCoul = 0.5∑n∑i∑j12piε0qiq j|r i j +n| , (2.4)where the first summation is over all periodic images which adds n = n1b1 +n2b2 + n3b3 to the position vector, r i j. n1,n2,and n3 are integers ranging from−∞ to +∞. Addition of n to r i j creates the periodic images of each set of par-ticles. Terms where i = j are ommitted when n = 0. Eq 2.4 decays at a rate of1/r, and is therefore only conditionally convergent. Truncation of eq 2.4 would20neglect a significant contribution to the potential energy of the simulation. A so-lution known as Ewald sums divides the Coulombic interations into short-rangeand long-range components.[98] The boundary conditions at infinity are conduct-ing. Short-range interactions between point charges are screened by a chargedistribution. The screening causes the interactions to decay faster, allowing fortruncation. Screening charges typically take the form of Gaussian distributionsof opposite sign to the charges they are screening. The resulting equation for theshort-range Coulombic interactions is[98]ushort =N∑i 6= jqiq jer f c(√αr i j)4piε0r i j, (2.5)where er f c is the complementary error function and α controls the width of theGaussian screening charge distributions. This solves the slow decay rate of eq2.4. However, the effect of the screening charges needs to be removed to pro-duce an accurate calculation of the total potential energy, creating a long-rangecontribution to the potential energy.In the Ewald sums, the Gaussian distributions are added back in, with oppositecharge, and the energy is calculated using a Fourier transform to produce a long-range contribution accounting for the truncation in eq 2.5.[98] The Ewald methodworks, but in most cases scales as the number of particles squared, N2, (Ewald canbe optimized to N3/2 in some cases) and becomes too slow to calculate in simula-tions larger than approximately 1000 particles. Faster methods of calculating thelong-range component have been developed.21GROMACS uses the smooth particle mesh Ewald method (SPME),[101–103]which scales as NlogN, to calculate the long-range portion of the Coulombic po-tential energy. This contribution to the potential is often referred to as the re-ciprocal energy because of the Fourier transforms used to calculate the energy inreciprocal space. SPME converts the distribution of charges to a discrete grid, ormesh, of charges and calculates the energy with a Fast Fourier transform method.Original charge positions are given via a structure factor,S(w) =N∑j=1q jexp(2piiw · s j), (2.6)where w = w1b∗1 + w2b∗1 + w3b∗1 are the reciprocal lattice vectors for the threeaxis of the simulation cell and w1,w2,and w3 are integers. The dot product ofw with the fractional coordinates s j, gives the coordinates of the particles in thesimulation and their images.To speed up the calculation, the structure factor is approximated by interpo-lating the charges onto mesh points.[101, 103] The interpolation is accomplishedwith a cardinal B-spline, Md(η). η are the fractional scaled coordinates of parti-cles in the simulation, and d is the degree interpolation. A cardinal B-spline is acontinuous spline with equally spaced knot points. The type of B-spline used isan Euler exponential spline which is differentiable, allowing for analytical calcu-lation of the energies and forces. The spline is defined by22M2(η) = 1−|η−1| for 0≤ η ≤ 20 for η < 0 or η > 2 , (2.7)where d = 2 is the lowest interpolation order, and a recursion formula is used ford > 2Md(η) =ηd−1Md−1(η)+d−ηd−1 Md−1(η−1). (2.8)Using the B-spline, the exponential in eq 2.6 is approximated byexp(2piiwiηi/Ki)≈ ci(wi)∞∑k=−∞Md(ηi− k)exp(2piiwik/Ki), (2.9)ci(wi) = exp(2pii(n−1)wi/Ki)[n−2∑k=0Md(k+1)exp(2piiwik/K1)]−1, (2.10)where K is the number of mesh points, and k is the location of the mesh points.Using the B-spline in eq 2.9 to assign charges to mesh points produces an approx-imate structure factorS˜(w) = c1(w1)c2(w2)c3(w3)F(Q)(w1,w2,w3), (2.11)where F(Q) is the Fourier transform of the charge array, Q. The charge arrayassigns charges, and the images of charges, to grid points, and is given by23Q(k1,k2,k3)=N∑i=1∑n1,n2,n3qiMd(η1i−k1−n1K1)Md(η2i−k2−n2K2)Md(η3i−k3−n3K3),(2.12)The reciprocal energy can then be approximated by[101, 103]u˜rec =12piV ∑w 6=0exp(−pi2w2/α2)w2S(w)S(−w) (2.13)=12piV ∑w 6=0exp(−pi2w2/α2)w2B(w1,w2,w3)F(Q)(w1,w2,w3)F(Q)(−w1,−w2,−w3)(2.14)whereB(w1,w2,w3) = |c1(w1)|2|c2(w2)|2|c3(w3)|2. (2.15)The above approximate reciprocal energy, u˜rec includes self interactions whichneed to be removed,[98]usel f =√α/piN∑i=1q2i4piε0, (2.16)so the full Coulombic potential energy is uCoul = ushort + u˜rec−usel f .242.3 Equations of MotionForces, F are calculated as the negative derivative of the energies, which allowsthe particles to be moved according to Newton’s equationsF i(r i j) =N∑j 6=i−∇u(r i j) = miai, (2.17)where F i(r i j) is the force on particle i, mi is the mass of the particle, ai is theacceleration, and ∇ is the gradient operator. Moving the particles is done by theLeap Frog algorithm which is derived from the Verlet algorithm,[98]r i(t +∆t) = 2r i(t)− r i(t−∆t)+∆t2 F i(t)mi , (2.18)where the position of the particle in the next time step, r i(t +∆t), is derived fromthe current position, the previous position, and the acceleration without requiringthe use of the particle’s velocity. To derive the Leap Frog algorithm, we beginwith the velocity, vi, of particle i half a time step before and after the current timestep,vi(t−∆t/2) = r i(t)− r i(t−∆t)∆t , (2.19)vi(t +∆t/2) =r i(t +∆t)− r i(t)∆t. (2.20)Updating the position of the particle is accomplished by rearranging eq 2.20 toobtain25r i(t +∆t) = r i(t)+∆tvi(t +∆t/2). (2.21)Updating the velocity is achieved from summation of eqs 2.19 and 2.20 followedby substitution of eq 2.18 to producevi(t +∆t/2) = vi(t−∆t/2)+∆t F i(t)mi . (2.22)The velocity is updated halfway between each time step in eq 2.22, and the posi-tion, r i, is updated every time step. In our simulations we use a 2 fs time step.Temperature coupling maintains the simulation at an average temperature, andis accomplished with the Nose´-Hoover thermostat.[104, 105] The Nose´-Hooverthermostat[104, 105] couples the simulation to a heat reservoir through a virtualparticle introduced into the Lagrangian of the simulation.[98]L =N∑i=1mi2s2v2i −u(rN)+Qm2s˙2− f kBT lns, (2.23)where the four terms on the right-hand-side give the kinetic and potential energiesof the particles in the simulation, and the kinetic and potential energies of theheat reservoir, respectively. mi, r i, and vi are the mass, position, and velocity ofparticle i in the simulation, respectively. f is the number of degrees of freedomin the simulation, kB is the Boltzmann constant, and T is the temperature of thesystem. s is the coordinate of the virtual particle, and Q the mass parameter ofthe virtual particle, which determines the coupling strength. This gives rise to athermodynamic friction coefficient, ζ , given by[105]26Qdζdt=∑imiv2i − ( f +1)kBT, (2.24)This thermodynamic friction coefficient changes the acceleration of the particlesin the system by[99]ai =F imi−ζ dr idt. (2.25)The result of the coupling is the simulation maintaining an average temperature.The Nose´-Hoover thermostat was used in all the simulations of this thesis,however, we did run a few tests using the velocity-rescaling thermostat[106] inkaolinite simulations, and we did not observe any significant differences in howice nucleated.Most simulations in this thesis use NV T conditions in which the number ofparticles, volume of the simulation cell, and the temperature remain constant.Occasionally, we allow the volume of the simulation cell to change, subject toconstant pressure conditions, NPT . To keep the pressure constant we use theParrinello-Rahman pressure coupling algorithm.[107, 108] Parrinello-Rahman pres-sure coupling is similar to the Nose´-Hoover thermostat. Both make changes to theLagrangian of the simulation and introduce new terms into the equations of mo-tion. In the Parrinello-Rahman barostat, the vectors of the simulation cell areintegrated according to the equations of motion allowing both the size and shapeof the cell to change. The simulation cell changes according to[99]27d2bdt2=VW−1b′−1(P−Pre f ), (2.26)where b′−1 is the inverse of the transpose of b, V is the volume of the cell givenby the determinant of b, W is a matrix determining the strength of the pressurecoupling, and P and Pre f are matrices of the current pressure and reference pres-sure, respectively. Pressure coupling has the following effect on the equations ofmotion for the particles in the simulation,ai =F imi−M dr idt, (2.27)M = b−1[bdb′dt+dbdtb′]b′−1. (2.28)Molecules, such as water, which are represented as a collection of sites withrigid internal structure, have their geometry preserved by the use of the LINCSconstraint algorithm.[109] In LINCS, at each time step, the constrained sites aremoved as if they were unconstrained and are then subject to two constraint stepsto return all the site-to-site distances (here referred to as “bonds”, even if no actualbond exists between the sites) in the molecule to their correct lengths. The firststep zeroes the projection of the new bonds onto the old bonds. The second stepshortens the bond length to what it should be. LINCS is performed on everybond, every time step.[99] In the case of water, constraints are applied on bothO-H bonds and the H-H distance to maintain the HOH angle.282.4 Force Fields and ModelsA number of water models have been created for molecular dynamics simulations.In this thesis we usually use two models to ensure our results are not stronglymodel dependent. All of the work includes the six-site water model[110] andeither the TIP4P/Ice[100] or the TIP4P/2005[72] model. Figure 2.1 shows thesix-site and TIP4P/Ice models. Both TIP4P/Ice and TIP4P/2005 share the samebasic design as the original TIP4P model,[111] modelling the oxygen atom ofwater with a LJ potential, placing positive charges on the hydrogen atoms, anda negative charge on what is called the “m-site”, which is located near the oxy-gen atom on the bisector of the HOH angle. Moving the negative charge off ofthe oxygen is required to qualitatively reproduce the phase diagram of water.[72]TIP4P/2005[72] is a reparametrization of the LJ parameters and partial chargesof the TIP4P model to create a model which accurately represents liquid water,and has a melting temperature of 250.5±3 K.[112] TIP4P/Ice[100] is a similarreparametrization, but designed to accurately model ice, and has a melting point of270±3 K.[112] Six-site water begins with the TIP4P model, and adds LJ terms tothe hydrogen atoms, and places small negative charges in the “lone pair” positions.The addition of “lone pairs” increases the propensity of the six-site model to posi-tion neighbouring water molecules more tetrahedrally compared to TIP4P models,i.e., making the angle formed by three water molecules closer to 109.5◦.[111] Themelting point of the six-site model is 289±3 K.[113] The m-site and lone pairs arerepresented as virtual particles in GROMACS.[99] Table 2.1 gives the parameters29for each of the water models.Table 2.1: Parameters for water models used in this thesis.Property six-site[110] TIP4P/2005[72] TIP4P/Ice[100]σO (nm) 0.3115 0.31589 0.31668εO (kJ/mol) 0.71485 0.7749 0.82165σH (nm) 0.0673 NA NAεH (kJ/mol) 0.115419 NA NAqH 0.477 0.5564 0.5897qm -0.866 -1.1128 -1.1794qLP -0.044 NA NAdOH(A˚) 0.98 0.9572 0.9572dOm(A˚) 0.23 0.1546 0.1577dOL(A˚) 0.8892 NA NAHOH6 (deg) 108 104.52 104.52LOL6 (deg) 111 NA NAIn this thesis water is simulated in contact with a variety of surfaces includingsilver iodide, kaolinite, feldspar, proteins, and gibbsite. Structures of these mate-rials are created using cif (crystallographic information files) files or pdb (proteindata bank) files obtained from various repositories.[114–116] Materials studiedinclude: silver iodide,[33, 114, 117] kaolinite,[118] feldspar,[119] gibbsite,[120]and a protein (PDB ID: 3ULT)[121]).Silver iodide interaction parameters were obtained from the original work ofHale and Kiefer.[35] Their model contained a polarization term, 1/r4, in additionto LJ and Coulombic terms. Polarization terms could not be implemented intoGROMACS because they are not pairwise additive, and were neglected when us-ing the Hale and Kiefer force field. We did perform some tests of the full Hale andKiefer force field, and we do not believe that excluding the polarization terms sig-30nificantly influences our results. Parameters used in this work are listed in Table2.2.Table 2.2: Parameters for the AgI and CLAFF[122] forcefields used in thisthesis. AgI values[35] are specific to the water model used.Atom Charge (e) σ (nm) ε (kJ/mol)Ag (TIP4P/Ice) 0.6 0.31732 5.93765I (TIP4P/Ice) -0.6 0.35132 7.67740Ag (six-site) 0.6 0.3225 7.32823I (six-site) -0.6 0.3565 9.47434Hydroxyl H 0.425 0.0 0.0Hydroxyl O -0.95 0.3166 0.650194Bridging O -1.05 0.3166 0.650194Substituted O -1.208 0.3166 0.650194Si 2.1 0.33020 7.70065e-06Al 1.575 0.4271 5.56388e-06K 1.0 0.3334 0.4184Calculations and a discussion of the polarization contribution to the silver-water and iodide-water interactions are given in Appendix B (see Figures B.1 andB.2). In Chapter 3, using infinite sheets of AgI, we tested the sensitivity of icenucleation to different interaction potentials by varying the partial charges on thesilver and iodide ions. We found that a wide range of potentials correspondingto partial charges ranging from 0 to ±0.8e nucleated ice efficiently. Figure B.2shows that when the induction term is included in the Hale–Kiefer force fieldthe potential continues to lie within (or very nearly so) the range of potentials forwhich we observe ice nucleation. Our observations indicate that the ice nucleatingability of the model AgI surfaces is not strongly dependent on details of the elec-trostatic water-AgI interactions. As discussed in detail in Chapter 3, the atomistic31geometry of the surface is the important factor.Kaolinite, feldspar, and gibbsite were all modelled using the CLAYFF forcefield.[122] This force field provides Lennard-Jones parameters and charges for anumber of ions which can be used to represent a large number of clays and miner-als. Parameters for the CLAFF force field were derived from DFT simulations forthe partial charges and the LJ terms were fitted to crystal structures. Validity of theforce field was tested by accurately producing several clay and mineral structures,such as kaolinite and pyrophyllite, and calculating power spectra which matchedexperimental infrared spectra.Several protein force fields have been developed,[123] and we used the am-ber03 force field.[124] which is a force field created for the simulation of proteinsand is already implemented in GROMACS. Parameters in the amber03 force fieldwere derived from DFT calculations and tested by simulating small proteins inwater.2.5 Simulation Setup2.5.1 Geometry of AgI, Kaolinite, Feldspar, and Gibbsite SlabSimulationsThe geometry of the simulation varied for each material, however, several projectsused a common two surface motif (Figure 2.2) to create a pair of infinite surfaces.Infinite slab simulations were used to determine the ice nucleation mechanism ofa particular surface of a material. One or more layers of a material were alignedin the X/Y plane of the simulation cell such that the periodic boundary conditions32create a continuous surface. Cleaving of crystal lattices can often produce sur-faces with a net dipole, which, combined with periodic boundary conditions, cancreate a long-range electrical field across the simulation cell capable of influenc-ing ice nucleation.[73] A second slab, mirroring the first in the Z direction,[88] isincluded in the simulation to cancel the spurious, long-range electric field. Watermolecules are placed between the slabs, and some empty space is left at the topand bottom of the simulation cell to separate the slabs because the simulation isalso periodic in the Z direction. Slab simulations of this form have been shown tobe effective in previous studies.[125] In our simulations we leave a total of 1 nmof vacuum space between the slabs (0.5 nm on each end of the simulation cell),which is the same as the short range cut-off used. Test simulations have been car-ried out with silver iodide where 2/3 of the simulation cell was vacuum,[126] andresults suggest that the difference in the amount of vacuum space does not havea significant effect on the results presented in this thesis. This configuration alsoprovides the ability to simulate two surfaces in one simulation. Particles formingeach surface are held fixed. Restricting particle movement in this manner in GRO-MACS means that the particle coordinates do not scale when the simulation cellchanges due to pressure coupling. This could lead to gaps in the surface if the boxexpands, or to the surface overlapping with its periodic images if the cell shrinks.Therefore, pressure coupling cannot be used with a mirroring slab configuration,and NV T conditions are employed.The number of water molecules between the slabs is adjusted to produce a bulkdensity between 0.92 and 0.96 g/ml, which has been found to be an effective range33WaterX/YZCrystal SlabVacuum(C)Figure 2.2: A diagram of the simulation cells used in much of this thesis.Opposing slabs of the material under consideration are shown in red,with the water between the slabs shown in blue. Vacuum areas at theend of the cell are white.for nucleating ice,[127] and gives room for ice to grow. We did carry out NPTsimulations of bulk liquid water at low temperatures (230 K for TIP4P/Ice and240 K for the six-site model), and in both cases the average bulk density at 1 barwas less than 0.95 g/mL, which is in the range we used. Therefore, although thelower densities used in NV T simulations might have some influence on nucleationrates (which we do not measure), we would not expect any significant influenceon the nucleation mechanism. Some NPT simulations were carried out with silver34iodide disks and plates as described in Chapter 4, and the results agreed with NV Tsimulations on the same systems.For each type of simulation, an initial run was carried out at 300 K. Coordi-nates were taken from these 300 K simulations every 0.5 ns to be used as initialcoordinates for production runs. Production runs were started at 300 K, equili-brated for 1 ns, then cooled to the desired temperature over 1 ns, and maintainedthere for the duration of the simulation. Temperatures of the simulations are ei-ther reported in Kelvin or in supercooling, ∆Ts, which is the number of degreesthe simulation is below the melting point of the water model used.2.6 AnalysisSimulations were analyzed using qualitative visual methods and quantitative meth-ods. Visual analysis was performed primarily with the VMD software,[128] andsometimes PyMOL[129] was also used. Most images were created with VMD.Density analysis was done for simulations with slab geometries using theGROMACS g density program. The density profiles show the average distribu-tion of atoms as a function of distance perpendicular to the surface. Densities arereported either normalized to the bulk concentration of the particular type of atomor in units of particles/nm3 to show hydrogen and oxygen atom densities on asimilar scale.Detection of ice in the simulations was done using the CHILL algorithm[130]implemented in a custom C++ code. The CHILL algorithm can detect both iceIh and ice Ic based on the arrangement of water oxygen atoms around a central35oxygen. As mentioned in Chapter 1, both forms of ice I are composed of hexag-onal rings, with Ic composed of only chair conformed rings. Ih is composed ofchair conformed layers stacked to produce boat conformed layers oriented per-pendicular to the chair conformed layers. CHILL classifies the bonds between awater oxygen atom and the four nearest water oxygen atoms as either staggered oreclipsed, based on the orientation of the four nearest water oxygen atoms aroundeach water molecule. As a result of the conformational differences, an oxygenin Ic will be surrounded by four other oxygen atoms in a staggered conformation,whereas an oxygen atom in Ih will by surrounded by three staggered oxygen atomsand one eclipsed oxygen atom. The CHILL algorithm detects this based on thethe following order parameter,qlm(i) = 0.254∑j=1Yl,m(r i j), (2.29)where qlm is a local order parameter around an oxygen atom, i, and its four nearestneighbours, j. Yl,m(r i j) is a spherical harmonic determined by the vector betweentwo neighbouring oxygen atoms, r i j. The l = 3 spherical harmonics are used,producing a set of seven complex values for each oxygen atom. Staggered andeclipsed configurations are detected by eq 2.30 which is a normalized dot productof vectors between oxygen atoms i and j,aq(i, j) =∑lm=−l qlm(i)q∗lm( j)(∑lm=−l qlm(i)q∗lm( j))1/2 (∑lm=−l qlm( j)q∗lm(i))1/2 , (2.30)36where q∗lm is the complex conjugate of qlm. A perfect staggered configuration willproduce a value of -1. An eclipsed configuration will produce a value of -0.11.In practice, a staggered configuration was defined as aq(i, j)<−0.75, an eclipsedbond was defined as having −0.25 < aq(i, j) < 0.1. Ice Ic would have four stag-gered bonds (one was allowed to have aq(i, j) < −0.7), and ice Ih would havethree staggered bonds and one eclipsed bond. Our values for aq(i, j) are a relaxeddefinition from the original CHILL algorithm.[130] The original algorithm wasdeveloped using the mW water model,[6] which has a greater propensity to adopttetrahedral configurations than the atomistic model used here.For the study of ice nucleation via AgI disks, it was necessary to determinethe size of the cluster of ice growing on a disk. One could not simply use the totalamount of ice in a simulation because the CHILL algorithm will detect some ice inthe solution due to random fluctuations, which has to be excluded from the clustersize determination. The DBScan algorithm[131] was employed to determine thesize of the cluster in the disk simulations. This algorithm detects clusters basedon the ability to connect particles by regions of some minimum density. Twoparameters are required by the algorithm: a radius, which was set to 3.8 A˚ here; aminimum cluster size, set to 4 in our calculations. The radius is based on oxygen-oxygen radial distribution functions for ice.[74] 3.8 A˚ captures the first peak in thedistribution representing the four closest water oxygen atoms.The DBScan algorithm, implemented in a C++ code, iterates through a list ofoxygen atoms classified as ice, which we obtain from the CHILL algorithm. If anatom is found to have four neighbouring ice oxygen atoms within the radius, then37a new cluster is created and the original atom and surrounding atoms are addedto the cluster. The cluster is then expanded by searching for atoms neighbouringthose in the cluster already.38Chapter 3A Molecular Mechanism of IceNucleation on Model AgI Surfaces“Science is magic that works.”— Kurt Vonnegut, Cat’s Cradle3.1 SummaryHeterogeneous ice nucleation at solid surfaces is important in many physical sys-tems including the Earths atmosphere. AgI is one of the best ice nucleating agentsknown; however, why AgI is such an effective ice nucleus is unclear. Usingmolecular dynamics simulations, we show that a good lattice match between iceand a AgI surface is insufficient to predict the ice nucleation ability of the sur-face. Seven faces modelled to represent surfaces of both β -AgI and γ-AgI, eachhaving a good lattice match with hexagonal and/or cubic ice, are considered, butice nucleation is observed for only three. Our model simulations clearly show39that the detailed atomistic structure of the surface is of crucial importance for icenucleation. For example, when AgI is cleaved along certain crystal planes twofaces result, one with silver ions and the other with iodide ions exposed as theoutermost layer. Both faces have identical lattice matches with ice, but in oursimulations ice nucleation occurred only at silver exposed surfaces. Moreover,although hexagonal ice is often the only polymorph of ice considered in discus-sions of heterogeneous ice nucleation, cubic ice was frequently observed in oursimulations. We demonstrate that one possible mechanism of ice nucleation byAgI consists of particular AgI surfaces imposing a structure in the adjacent waterlayer that closely resembles a layer that exists in bulk ice (hexagonal or cubic).Ice nucleates at these surfaces and grows almost layer-by-layer into the bulk.3.2 MethodsSeven surfaces of AgI are tested for their ice nucleating ability: the silver andiodide exposed β -(0001), γ-(001), and γ-(111), and the β -(101¯0), which has nosilver or iodide exposed forms, and thus we refer to it as Ag+I exposed. Simu-lations involved the six-site[110] and the TIP4P/Ice[100] models of water. Mostsix-site simulations were run at 270 K, and TIP4P/Ice simulations at 250 K, soboth were supercooled by about 20 K. Dimensions of the simulations used forAgI slabs are given in Table 3.1. Lists of simulations used in this chapter are pro-vided in Tables A.1, A.2, A.3, and A.4, and a discussion of the effect of mirroringslabs is provided in Appendix A and in Figures 3.1 and 3.2. Results presented inthis chapter are for the six-site model. Results for the TIP4P/Ice model are similar40to the six-site model and are provided in Appendix A.Table 3.1: Summary of the surfaces considered. XY dimensions, the num-bers of water molecules, and bulk liquid densities of liquid water aregiven. All surfaces except β -AgI(101¯0) were simulated with both sil-ver and iodide ions exposed, and the surfaces are 5 nm apart measuredfrom the closest ions.Surface X/Y Dimensions (nm) Number of water molecules Bulk densities (g/mL)six-site TIP4P/Ice six-site TIP4P/Iceβ -0001 5.568/5.51 4780 4610 0.966 0.933β -101¯0 6.008/5.51 5000 5000 0.926 0.932γ-111 5.515/5.572 4780 4610 0.964 0.928γ-001 5.849/5.849 5300 5300 0.931 0.9343.3 PolarizationSince the model AgI lattice has point charges on the Ag and I sites, it is interestingto calculate polarization profiles of the ice that is formed. Significant polarizationmight indicate that ice forms under the influence of an electric field created bythe surface, and thus nucleation may not be a result of the silver iodide surfaceserving as a template for ice. An average dipole order parameter or polarization isgiven by:[75]〈m(z)〉=〈1N(z)µwN(z)∑i=1µ i〉(3.1)where µ i is the dipole vector of water molecule i, µw is the dipole moment of thewater model considered, z is the distance from the surface, and N(z) is the numberof water molecules in a rectangular volume element of width 0.001 nm centred atz. The components of 〈m(z)〉 for the three surfaces that nucleated ice are shown in41<m(z)>-0.500.511.522.53z (nm)2 3 4 5 6ABC1Figure 3.1: Polarization profiles for the six-site model. Panels A, B, and Cshow results for Ag exposed β -AgI(0001), Ag exposed γ-AgI(111),and Ag exposed γ-AgI(001), respectively. The red, green, and bluecurves are the X, Y and Z components, respectively. Note the shiftedscale in Panels A and B.Fig. 3.1 (six-site model) and Fig. 3.2 (TIP4P/Ice). These profiles were obtainedfor frozen samples. Note that the Z component (blue) reverses sign because thesurfaces face in opposite directions. We see that for both models the polarization isonly significantly nonzero in the Z direction, and that it is very weak except nearthe surface. The polarization oscillations are a little stronger for the TIP4P/Icemodel (Fig. 3.2) likely because its dipole moment (2.426 D[100]) is larger thanthat of the six-site model (1.89 D[110]). However, based on earlier work,[75] forboth water models the surface-induced polarization is much too weak to cause icenucleation.42<m(z)>-0.500.511.522.53z (nm)2 3 4 5 6ABC1Figure 3.2: Polarization profiles for the TIP4P/Ice model. Details are as inFig. 3.1.3.4 ResultsSeven different silver iodide surfaces were considered (Table 3.2). The differentsurfaces were generated by cleaving either hexagonal AgI (referred to as β -silveriodide) or face centred cubic (FCC) AgI (referred to as γ-silver iodide). Bothforms of AgI can exist under atmospheric conditions. If we cleave a β -AgI(0001)plane, a γ-AgI(111) plane, or a γ-AgI(001) plane, two faces are produced, onewith sliver ions exposed and one with iodide ions exposed. We refer to these twofaces as silver exposed and iodide exposed, respectively. Cleaving a β -AgI(101¯0)plane produces two identical faces, so this face does not have silver or iodide ex-posed versions. We refer to this surface as the silver+iodide exposed β -AgI(101¯0)face.All surfaces considered have a small mismatch with hexagonal ice and/or cu-43Table 3.2: Summary of the AgI faces considered, the associated mismatcheswith Ih and Ic, and the polymorph(s) of ice formedAgI polymorph AgI surface Mismatch (%) Ice formedIh IcβAg exposed (0001) 1.55 2.00 Ih+IcI exposed (0001) 1.55 2.00 noneAg+I exposed (101¯0) 2.03 13.45 noneγAg exposed (111) 1.77 2.22 Ih+IcI exposed (111) 1.77 2.22 noneAg exposed (001) 41.61 2.04 IcI exposed (001) 41.61 2.04 nonebic ice (Table 3.2). The lattice mismatch can be quantified with the equation[1]δ =aAgI−aiceaice×100% , (3.2)where δ is the lattice mismatch, and aice and aAgI are the lattice constants of iceand the silver iodide surfaces, respectively. Although previous researchers oftenonly considered the lattice mismatch between hexagonal ice and AgI surfaces,[1,29, 132, 133] here we consider the lattice mismatch with both hexagonal ice andcubic ice, for reasons that will become clear below.Even though all AgI surfaces considered have a small lattice mismatch with Ihand/or Ic, ice formation was only observed for three of the surfaces (Table 3.2, lastcolumn). For the three surfaces that nucleated ice, ice nucleated and grew nearlyas soon as the liquid was supercooled by 20 K (Figures 3.3 & A.1).Shown in Figure 3.4 are snapshots from simulations for two surfaces that suc-cessfully nucleated ice at a supercooling of 20 K. Snapshots for the third surface44ice01,0002,0003,0004,0005,000time at 270K (ns)0 20 40 60 801Figure 3.3: The number of ice molecules as a function of time detected bythe CHILL algorithm in simulations on three AgI surfaces: Ag ex-posed β -AgI(0001) (red), Ag exposed γ-AgI(111) (green), and Ag ex-posed γ-AgI(001) (blue). The solid lines and dotted lines represent thenumber of Ic and Ih molecules, respectively.that successfully nucleated ice under the same conditions are shown in Figure3.5. The CHILL algorithm[130] was used to determine whether a particular watermolecule was part of the liquid, Ih, or Ic phase. Note that this algorithm does notrecognize water molecules at the silver iodide surface as ice, so they retain the liq-uid colour code, however, these water molecules form part of the initial ice layerafter nucleation has occurred.Figure 3.4A-C shows snapshots for the silver exposed β -AgI(0001) face. Forthis face both Ih(0001) and Ic(111) planes grow parallel to the silver iodide sur-face. As remarked earlier, the Ih(0001) and Ic(111) faces are almost identical andcan connect with little lattice strain.[23] This stacking disorder has been noted inprevious work.[24, 25]45Figure 3.4: Snapshots of ice growth at different times. Panels A-C are forthe silver exposed β -AgI(0001) face at 0, 11 and 36 ns, respectively.Panels D-F are for the silver exposed γ-AgI(001) face at 0, 3 and 6ns, respectively. Only the first layer of each AgI slab is shown. Silverand iodide ions are silver and green, respectively. Oxygen atoms ofliquid water, Ic, and Ih are red, yellow and blue, respectively, and allhydrogen atoms are black. Note that panels D-F have been rotated toshow the (110) plane of Ic, and this is why the water density can appearnonuniform in those snapshots.46A B CFigure 3.5: Snapshots of ice growing on the Ag exposed γ-AgI(111) surface.Panels A-C are taken at 0 ns, 7 ns, and 15 ns, respectively. Silverand iodide ions are silver and green, respectively. Oxygen atoms ofliquid water, Ic, and Ih are red, yellow and blue, respectively, and allhydrogen atoms are black. The panels are rotated as in Fig. 1D-E.In Figure 3.4A-C, at ∼ 11 ns the ice layers growing from both silver iodidesurfaces meet, the upper crystal then continues to grow at the expense of the lowerone, in what appears to be an Ostwald ripening process. In the simulation, the icegrowth eventually stops with some liquid water remaining in the simulation cell.This is because as the ice grows the liquid is compressed until no further icegrowth is possible.3.4.1 Why do some model AgI surfaces nucleate ice andothers not?Comparison of simulations for Ag exposed and I exposed β -AgI(0001)To better understand why some AgI surfaces nucleate ice and others do not, eventhough they all have a small lattice mismatch with at least one ice polymorph,we compare the microscopic details of the water-surface interaction for the sil-ver exposed β -AgI(0001) and the iodide exposed β -AgI(0001) faces (Figure 3.6).47Both faces have identical lattice mismatches but only the silver exposed face nu-cleates ice. On both the silver and iodide exposed β -AgI(0001) faces, water formshexagonal rings (Figure 3.6A & B). This finding is consistent with previous stud-ies that have shown that the β -AgI(0001) surface is effective at organizing watermolecules on its surface into hexagonal rings.[34–43]For the silver exposed surface (Figure 3.6A & C), the oxygen atoms of the wa-ter molecules have a chair conformation, which resembles the chair conformationin the ice bilayer of Ih(0001) and Ic(111) (compare panels C and E in Figure 3.6).It is apparently this striking conformational similarity that makes this surface anexcellent ice nucleus. The chair conformation is caused by water molecules oc-cupying two different sites at the surface: directly above iodide ions lower in theAgI bilayer, or over interstitial sites between both ions. Water molecules occupy-ing interstitial sites sit closer to the surface than water molecules over iodide ionsand thus form the chair conformation.In contrast, for the iodide exposed surface (Figure 3.6B & D) the hexago-nal rings form a more coplanar structure. On the iodide exposed face, all watermolecules sit over equivalent lattice sites, such that all oxygen atoms are largelycoplanar. This planar layer of hexagonal rings does not match any surface of ice,which is likely why we observe no ice nucleation on this face.The difference between these two faces can also be seen in the normalizeddensity profiles of oxygen and hydrogen, as functions of distance from the sur-face (Figure 3.7, Figure A.2 for TIP4P/Ice). For the silver exposed face densityprofiles are given for both ice (Figure 3.7A) and liquid water (Figure 3.7B) near48A BC DEFigure 3.6: Snapshots of liquid water on β -AgI(0001) surfaces. Panel A is atop view of water (supercooled by 4 K) on the silver exposed face, andPanel C is a side view of a single hexagonal ring on this surface. PanelsB and D are corresponding views of water (supercooled by 20 K) onthe iodide exposed face. Panel E is a side view of a chair structure fromice. Silver ions, iodide ions, oxygen atoms, and hydrogen atoms aresilver, green, red and black, respectively. The yellow lines highlightthe hexagonal arrangement of the water on the surface.the surface. For ice (Figure 3.7A), the oxygen profile shows a repeating doubletindicative of the formation of bilayers with hexagonal rings in the chair conforma-tion. The hydrogen density profile for ice consists of repeating triplets, the centralpeak due to hydrogen atoms that are hydrogen bonded with water molecules in49 ρ/ρ o02468101214distance from surface (nm)0 0.5 1ABC1Figure 3.7: Normalized density profiles of oxygen atoms (red) and hydro-gen atoms (black) of water as functions of distance from β -AgI(0001)surfaces. Panels A and B are for ice (supercooled by 20 K) and liquidwater (supercooled by 4 K), respectively, near the silver exposed face.Panel C is for liquid water (supercooled by 20 K) near the iodide ex-posed face. On the distance scale, zero (vertical green line) coincideswith the position of the centre of the outermost ion of the face. ρ0 isthe average density of oxygen or hydrogen atoms far from the surface.Note the shifted scale in Panels A and B.the same bilayer. The less intense peaks on either side indicate hydrogen bondingto ice layers above and below the bilayer, or in the case of the initial peak, to thesilver iodide surface. Liquid simulations (Figure 3.7B) show that near the surfacethe structure of water is very similar to that of ice. It is apparent that the silverexposed β -AgI(0001) surface imposes an icelike structure in the adjacent waterlayer. This icelike surface layer then allows the next layer of ice to form and soon as ice continues to grow from the surface. In all likelihood, this is why thisparticular AgI face is such an effective ice nucleus.The iodide exposed face of β -AgI(0001) did not nucleate ice. The density50profiles shown for this face are for the liquid state (Figure 3.7C), and are to becompared with those for liquid near the silver exposed face (Figure 3.7B). It canbe seen that the water structure near the iodide exposed face differs significantlyfrom the silver exposed case. On the iodide exposed face the density profiles forboth the oxygen and hydrogen atoms exhibit single rather than split peaks. Thehydrogen peak is nearly coincident with that of oxygen indicating that the hy-drogen and oxygen atoms are essentially coplanar. The water density profiles forthe iodide exposed surface suggest a single planar layer of water molecules ratherthan a bilayer, or, in other words, the hexagonal rings at the surface (Figure 3.6B)are planar rather than in chair conformations. A similar density profile, with alarge single surface peak, was also observed in an earlier simulation of water onthis iodide exposed face.[40] The failure of the iodide exposed face to order waterinto an icelike bilayer likely undermines its ability to nucleate ice.The results discussed above were all obtained with charges of±0.6e on the sil-ver and iodide sites; we shall refer to this case as the “normally” charged surface.However, as noted above, simulations were performed for a range of charges, andthese give additional insight into how the silver exposed face organizes water intochair configurations. For the silver exposed β -AgI(0001) face, ice nucleation wasobserved for charges ranging from zero to ±0.8e, but not for ±1.0e (Table A.4).This behaviour can be understood by examining the charge dependence of thedensity profiles of liquid water near the silver exposed β -AgI(0001) face; profilesfor charges of zero and±1.0e are shown in Figure 3.8 (Figure A.6 for TIP4P/Ice).We note that the profiles for the zero charge case (Figure 3.8A) resemble those51 ρ/ρ o012345678distance from surface (nm)0 0.5 1AB1Figure 3.8: Normalized density profiles of oxygen atoms (red) and hydrogenatoms (black) of liquid water as functions of distance from the Agexposed β -AgI(0001) surface. Panel A is the zero charge case (watersupercooled by 4 K). Panel B is for AgI carrying charges of ±1.0e(water supercooled by 20 K). The green vertical line at zero coincideswith the position of the centre of the outermost ions, and ρ0 is theaverage density of oxygen or hydrogen atoms far from the surface.Note the shifted scale for Panel A.of liquid water near the normally charged face (Figure 3.7B). Thus, the orderingof water at this surface appears to be mainly due to the LJ interactions, which areall that remain in the zero charge limit. On the other hand, for charges of ±1.0e,the density profiles (Figure 3.8B) show significant deviations from those of thenormally charged surface. In particular, the high sharp peaks suggest that watermolecules at this surface are more tightly held, with the oxygen atoms occupy-ing a more planar arrangement. This distortion of the surface layer inhibits icenucleation.The γ-AgI(111) faces are very similar to the corresponding faces of β -AgI(0001)52A BC DFigure 3.9: Snapshots of water on the γ-AgI(111) surfaces. Panel A is atop view of water (supercooled by 4 K) on the silver exposed face,and Panel C is a side view of a single hexagonal ring on this surface.Panels B and D are corresponding views of water (supercooled by 20K) on the iodide exposed face. Colours for silver ions, iodide ions,oxygen atoms and hydrogen atoms are silver, green, red, and black,respectively. Yellow lines highlight the hexagonal arrangement of thewater on the surfaceand order water in an near identical manner (Figures 3.9 & 3.10 & 3.5 & A.4).Thus, the discussion above for β -AgI(0001) applies equally to the γ-AgI(111)case.53 ρ/ρ o02468101214distance from surface (nm)0 0.5 1ABC1Figure 3.10: Density profiles of oxygen (red) and hydrogen (black) atomsfor the six-site model near the γ-AgI(111) surface. Panels A and Bare for ice (supercooled by 20 K) and liquid water (supercooled by4 K), respectively, near the silver exposed face. Panel C is for liquidwater (supercooled by 20 K) near the iodide exposed face. On thedistance scale, zero (vertical green line) coincides with the positionof the centre of the outermost ion of the surface. ρo is the averagedensity far from the surface. Note the shifted scale in Panels A andB.Comparison of simulations for Ag exposed and I exposed γ-AgI(001)It is also interesting to compare the microscopic details of the water-surface in-teractions for the silver exposed face of γ-AgI(001) with the iodide exposed faceof γ-AgI(001). Both surfaces have identical lattice mismatches to Ic, but only thesilver exposed face led to ice nucleation. Configurational snapshots of the initialwater layer on the silver and iodide exposed faces of γ-AgI(001) are shown in Fig-ure 3.11A & B. These snapshots are for liquid water in contact with the surface.We note that for both faces the oxygen atoms occupy positions consistent with a54cubic ice lattice. However, the molecular orientations are different; most notably,for the silver exposed surface, protons of the water molecules are clearly directedinward towards the surface, similar to the structure of water in an Ic lattice (Fig-ure 3.11E). This fine detail appears to be important for the nucleation of ice. Incontrast, for the iodide exposed surface, all of the hydrogen atoms are directedoutward away from the surface, which does not match any surface of ice. Thismismatch in the fine structure is likely why we observe no ice nucleation on thisface.The differences between the two faces can also be seen in the density profilesshown in Figure 3.12 (Figure A.3 for TIP4P/Ice). For the silver exposed face den-sity profiles are given for both ice (Figure 3.12A) and liquid water (Figure 3.12B).For ice (Figure 3.12A), the oxygen profile presents a series of evenly spaced peaksthat are from the layers of the cubic lattice of Ic. The peaks in the hydrogen atomprofile do not coincide with those of oxygen, indicating that the hydrogen atomsare bonding to water molecules in other layers or to the silver iodide surface, butnot to water molecules in the same layer. Profiles for liquid water near the silverexposed face (Figure 3.12B) clearly show that icelike structure is imprinted on thewater at the surface, just as it is for silver exposed β -AgI(0001) case, allowing thenucleation and growth of subsequent ice layers.Density profiles for the iodide exposed face do not match the ice structureas well as the silver exposed case (Figure 3.12C). This is particularly true of thehydrogen atom profile, which clearly shows the absence of protons closer to thesurface than oxygen atoms. Thus, just as for β -AgI(0001), the surface induced55A BC DEFigure 3.11: Snapshots of liquid water on γ-AgI(001) surfaces. Panel A isa top view of water (11 K above the melting point) on the silver ex-posed face, and Panel C is a tilted view of a single FCC-like arrange-ment on this surface. Panels B and D are corresponding views ofwater (supercooled by 20 K) on the iodide exposed face. Panel E isa tilted view of a FCC arrangement in Ic. Silver ions, iodide ions,oxygen atoms, and hydrogen atoms are silver, green, red and black,respectively. The yellow lines highlight the FCC-like arrangement ofthe water on the surface.56 ρ/ρ o02468101214distance from surface (nm)0 0.2 0.4 0.6ABC1Figure 3.12: Normalized density profiles of oxygen atoms (red) and hydro-gen atoms (black) of water as functions of distance from γ-AgI(001))surfaces. Panels A and B are for ice (supercooled by 20 K) and liq-uid water (6 K above the melting point), respectively, near the silverexposed face. Panel C is for liquid water (supercooled by 20 K) nearthe iodide exposed face. On the distance scale, zero (vertical greenline) coincides with the position of the centre of the outermost ion ofthe face. ρ0 is the average density of oxygen or hydrogen atoms farfrom the surface. Note the shifted scale in Panels A and B.organization of the initial layer or two of water molecules is crucial for ice nucle-ation.Some simulations with different surface charges were also done for the silverexposed γ-AgI(001) face (Table A.4). Ice nucleation was observed for chargesranging from ±0.2e to ±0.6e, but not for zero and ±0.8e. Density profilesobtained for zero charge and ±0.8e are plotted in Figure 3.13 (Figure A.7 forTIP4P/Ice). Both profiles differ considerably from the ±0.6e result for this sur-face (Figure 3.12B). At zero charge (Figure 3.13A), the water molecules appearto have moved away from the surface to some extent, indicating that in this case57 ρ/ρ o01234567distance from surface (nm)0 0.5 1AB1Figure 3.13: Normalized density profiles of oxygen atoms (red) and hydro-gen atoms (black) of liquid water (supercooled by 20 K) as functionsof distance from the Ag exposed γ-AgI(001) surface. Panel A is thezero charge case. Panel B is for AgI carrying charges of ±0.8e. Thegreen vertical line at zero coincides with the position of the centre ofthe outermost ions, and ρ0 is the average density of oxygen or hy-drogen atoms far from the surface. Note the shifted scale for PanelA.the Coulombic interactions do significantly influence how water is organized atthe surface. Conversely, the profiles for ±0.8e (Figure 3.13B) show that wateris more tightly held by the higher charge. Additionally, at the surface the watermolecules are now oriented such that both protons are directed inward towardsiodide ions (Figure 3.14), which likely inhibits ice nucleation.Simulations for the Ag+I exposed β -AgI(101¯0) surfaceAs discussed above, when β -AgI(101¯0) is cleaved only one type of surface isformed and we refer to this as Ag+I exposed, since both silver and iodide ions arein the plane of the surface. Despite the good lattice match to Ih, ice nucleation58Figure 3.14: Snapshot of TIP4P/Ice water near the Ag exposed γ-AgI(001)face with charges of ± 1.0e on the AgI (supercooled by 20 K).Colours for silver ions, iodide ions, oxygen atoms, and hydrogenatoms are silver, green, red, and black, respectively.was not observed for this particular silver iodide face. Configurational snapshotsshow that the first layer of water on this face does not possess any significanticelike order (Figure 3.15) and density profiles do not show a strong icelike pattern(Figures 3.16 & A.5). Thus, it appears that the atomistic details of this surface donot favour ice nucleation and growth.59Figure 3.15: The first layer of water on the Ag+I exposed β -AgI(101¯0) face(supercooled by 20 K). Colours for silver ions, iodide ions, oxygenatoms, and hydrogen atoms are silver, green, red, and black, respec-tively. The yellow and light blue lines highlight the two orientationsof the boat conformations made by the silver iodide surface; yellowis opening outward and blue is opening inward towards the surface.60 ρ/ρ o00.511.522.53distance from surface (nm)0 0.5 1 1.51Figure 3.16: Density profiles of oxygen (red) and hydrogen (black) atomsfor six-site water near the Ag+I exposed β -AgI(101¯0) surface (su-percooled by 20 K). On the distance scale, zero (vertical green line)coincides with the position of the centre of the outermost ion of thesurface. ρo is the average density far from the surface.613.5 ConclusionsBased on our model simulations and the detailed analysis discussed above, a pos-sible mechanism for ice nucleation by silver iodide consists of some silver iodidesurfaces imposing a structure in the adjacent water layer that closely resembles alayer that exists in bulk ice (either Ih or Ic). Ice nucleates at these surfaces andgrows outward forming successive layers into the bulk. This mechanism does notrule out defects being important for ice nucleation by silver iodide, but does indi-cate that defects are not necessary in order to explain why silver iodide particlescan be highly effective ice nuclei. However, as noted by Fraux and Doye,[83] onemust be cautious in extrapolating from model simulations to real AgI surfaces. AllAgI faces for which we observed ice nucleation have surface dipole moments, andhence some surface reconstruction is expected.[134] Clearly, such effects, as wellas any other possible influences of thermal motion, are not included in our rigidlattice model of AgI. Finally, we speculate that a templating mechanism similarto that observed for our model AgI surfaces might apply for other effective icenuclei, such as aliphatic alcohol monolayers, and proteins from the bacteria Pseu-domonas syringae.62Chapter 4Simulations of Ice Nucleation byModel AgI Disks and Plates“People have to talk about something just to keep their voice boxes inworking order, so they’ll have good voice boxes in case there’s everanything really meaningful to say.”— Kurt Vonnegut, Cat’s Cradle4.1 SummarySilver iodide is one of the most effective ice nuclei known. We use moleculardynamics simulations to investigate ice nucleation by AgI disks and plates withradii ranging from 1.15 to 2.99 nm. It is shown that disks and plates in this sizerange are effective ice nuclei, nucleating bulk ice at temperatures as warm as 14K below the equilibrium freezing temperature, on simulation time scales (up to afew hundred nanoseconds). Ice nucleated on the Ag exposed surface of AgI disksand plates. Shortly after supercooling an ice cluster forms on the AgI surface. The63AgI-stabilized ice cluster fluctuates in size as time progresses, but, once formed,it is constantly present. Eventually, depending on the disk or plate size and the de-gree of supercooling, a cluster fluctuation achieves critical size, and ice nucleatesand rapidly grows to fill the simulation cell. Larger AgI disks and plates supportlarger ice clusters and hence can nucleate ice at warmer temperatures. This workmay be useful for understanding the mechanism of ice nucleation on nanoparticlesand active sites of larger atmospheric particles.4.2 MethodsFinite-size disks and plates of different shape are carved from four layer thicksheets of β -AgI,[33] such that a (0001) face forms the top (silver exposed) andbottom (iodide exposed) of each disk or plate. Top views of the disks and platesconsidered are shown in Figure 4.1. Note that all disks and plates are four layersthick and are electrically neutral containing equal numbers of Ag and I ions. Disksand plates labelled A-E were obtained by retaining all atoms within a specifiedcircumscribing radius, which naturally leads to hexagonal shapes for the smallerplates (A-C) and to roughly circular shapes for the larger disks (D,E), as can beseen in Figure 4.1. The rectangular plate (labelled Dr) has the same number ofatoms as disk D and is considered to probe for possible effects of shape on thenucleation properties. The circumscribing radii of the disks and plates togetherwith the rectangular dimensions of plate Dr are given in Table 4.1.In all simulations, with one exception (discussed later), the simulation cellcontained a single AgI disk or plate, centred in the XY plane with the (0001) face64Figure 4.1: Top view of the AgI disks and plates considered. Silver andgreen spheres represent silver and iodide ions, respectively.65Table 4.1: Summary of sytems considered. System C2 contains 2 mirroringplates, and system Dr is a rectangular plate. X/Y/Z gives the dimensionsof the simulation cell. ∆Ts reports the supercooling required to nucleateice. Densities are for the six-site model at 300 K.System Disk radius (nm) Water molecules X/Y/Z (nm) ∆Ts (K) Density (g/mL)A 1.15 5350 5/5/7 35 0.947B 1.61 7620 6/6/7 26 0.947C 2.07 7450 6/6/7 22 0.952C2 2.07 17650 6/6/16.44 21 0.947D 2.53 11870 7/7/8.222 16 0.945Dr 4.507×4.229 11870 7/7/8.22 16 0.945E 2.99 17150 8/8/9 14 0.953perpendicular to the Z axis. Periodic boundary conditions are applied in all direc-tions so the actual location of the disk or plate along the Z axis has no influenceon the physical behaviour of the system; however, for purposes of analysis anddepiction the disk or plate is located near the bottom of the cell with the silver ex-posed surface on which ice nucleates directed upward in the positive Z direction.The cell dimensions and the number of water molecules included for the differentsystems are given in Table 4.1.Most simulations were carried out under NV T conditions with water densitieslying in the range 0.94 to 0.96 g/mL (Table 4.1); however, to check that our re-sults are not strongly influenced by the relatively low densities (closer to ice thanto liquid water) or the constant volume condition, simulations for systems A, B,C, and D (Table S2) were carried out under NPT conditions with the pressurecontrolled at 1 bar using the Parrinello-Rahman barostat.[107, 108] In these NPTsimulations, ice nucleated on the surfaces via the same mechanism and in simu-66lation times comparable to the NV T results (see Appendix B). Thus, at least forthe conditions we consider, the simulation results are not sensitive to the ensem-ble employed. A list of all the simulations employed in Chapter 4 is provided inTables B.1, B.2, and B.3.The configuration labelled C2 contained two mirrored plates located 10 nmapart in the simulation cell. This system was used to check that spurious electricfields were not influencing our single disk or plate simulations, as can be the casefor slabs of charged particles that are periodically infinite in two dimensions, butof finite thickness.[88] We also carried out simulations (for systems C and D)with no charges on the sliver and iodide ions. Ice nucleation occurred in thesesimulations, just as it does when the full force field (with charges of ±0.6e) isused (see Appendix B). This provides a further demonstration that surface fieldsare not strongly influencing our results.4.3 PolarizationIn our simulations, we might expect the AgI disks and plates to give rise toan electric field. Large electric fields that highly polarize water can induce icenucleation.[75] Therefore, to check this possibility we calculated polarizationprofiles[75]〈m(z)〉=〈1N(z)µwN(z)∑i=1µ i〉, (4.1)67where 〈m(z)〉 is the net polarization (mean reduced dipole moment) in a rectangu-lar volume element of width 0.01 nm at z, N(z) is the number of water moleculesin the volume element, µw is the molecular dipole moment of the water model,and µ i is the dipole vector of molecule i. The components of 〈m(z)〉 are plottedin Figure 4.2.The plots show that the net polarization of water is too small to significantlyinfluence ice nucleation. This is confirmed by the fact that systems C and C2 (inC2 the fields from the disks cancel, such that there is no net polarization awayfrom the disk surfaces) give very similar nucleation results.4.4 RationalIn this chapter, we investigate ice nucleation on AgI disks and plates with nanome-ter dimensions. This investigation provides insight into the mechanism of ice nu-cleation on nanoparticles. In addition our simulations allow us to directly explorethe size, nature, and relevance of critical clusters during heterogeneous nucleation.We note that such free-standing AgI disks and plates may not exist in real physicalsituations; however, it may be possible for a larger silver iodide particle to havesmall patches that are of similar size and structure to those studied here. Hence,our studies may provide insight into the ice nucleation properties of larger sil-ver iodide particles. More generally, atmospheric particles, which typically rangefrom 10 nm to 10 µm in size,[3] may present many different crystal faces. As a re-sult the surface of atmospheric particles may also contain crystalline patches withnanometer dimensions similar to those studied here. For convenience we refer to68Figure 4.2: Polarization profiles of water at ∆Ts values 1 K above where icenucleation was observed. Red, green, and blue represent the X , Y andZ components, respectively. The system label and the ∆Ts value aregiven in the legend. In all systems except C2, the silver exposed faceof the AgI disk is directed outward along the positive Z axis. For C2,the silver exposed faces are directed towards each other, such that anylong-range fields cancel.69these patches as active sites. Our studies may also be useful for understanding themechanism of ice nucleation on active sites of this type.4.5 Results and DiscussionWe observe that ice nucleation on disks and plates with nanometer dimensionsoccurs in two stages. Initially, several layers of ice form on the disk or plate sur-face. These layers fluctuate growing and shrinking as the simulation proceeds, butthe surface ice cluster never melts, and the number of permanently frozen watermolecules can represent a significant fraction of the required critical nucleus. Fora particular disk or plate, if the temperature is low enough, an ice cluster fluctua-tion eventually matches the critical size, and bulk ice nucleates and rapidly growsto fill the simulation cell. As the temperature is increased, larger disks are neededto create the larger critical clusters that are required to achieve nucleation. Wher-ever it is possible and meaningful, we compare our estimates of critical clustersize with theory, previous simulations, and experiment.In our simulations ice clusters form only on the silver exposed face of the AgIdisks and plates, consistent with previous work showing that this face nucleatesice very effectively.[83, 135] If the surface-induced ice cluster achieves a criticalsize (dependent on ∆Ts), then ice nucleates and grows throughout the simulationcell.Several configurational snapshots of an ice cluster on the silver exposed sur-face of disk D nucleating bulk ice at ∆Ts = 16 K are shown in Figure 4.3. We notethat at earlier times in the simulation (30-91 ns in Figure 4.3) the ice cluster is con-70Figure 4.3: Series of snapshots of an ice cluster nucleating bulk ice on AgIdisk D. At each time, side and top views of the ice cluster are shown.Silver and iodide ions are shown in silver and green, respectively. Inthe first and third columns, black, blue, and orange represent hydrogenatoms, the oxygen atoms of Ih, and the oxygen atoms of Ic, respec-tively. In the second and fourth columns, yellow and purple representthe oxygen atoms of ice inside and outside the cylindrical volume de-fined by the disk radius, respectively. Liquid water molecules are notshown.tained within the cylindrical boundary of the disk and is roughly hemispherical inshape. At 91 ns the cluster has achieved critical size and begins to extend beyondthe cylindrical disk boundary (see later for further discussion on this point), and itthen grows primarily horizontally eventually filling the simulation cell. As iden-tified by the CHILL algorithm, the ice on this disk (and the other disks and plates71considered) is a mixture of Ih and Ic, or Isd .[26] The c–axis of Ih is perpendicularto the disk surface. Thus, the initial ice growth is on the basal face as the clustergrows taller. Later, horizontal growth occurs on what would be prism faces if theice was completely Ih.For each system considered, the smallest ∆Ts for which bulk ice nucleationwas observed is given in Table 4.1. For the disks and plates (Figure 4.1AD) ∆Ts isplotted versus the disk or plate radius in Figure 4.4 and compared with previousresults. The black horizontal line in Figure 4.4 is the minimum observed ∆Ts (8K) required to nucleate ice on infinite silver iodide slabs (Chapter 3). The trendin our results indicates that the infinite surface limit is being approached as thedisk or plate radius increases. We note (Table 4.1) that the rectangular surface(Dr) nucleated ice at the same ∆Ts as the disk D, suggesting that the particularshape of the exposed AgI surface is not a major factor influencing ice nucleation.Also, the mirrored plate system C2 nucleated ice only one degree warmer than thesingle disk or plate system, indicating that the small electric fields and resultingpolarization of the water produced by single disks or plates (see Figure 4.2) donot significantly influence the nucleation temperature. It is possible that the onedegree difference found for the C2 case is due to the ice-like clusters on bothplates interacting, making it easier to achieve critical cluster size. Additionally,note that the C and D surfaces nucleate ice at the same temperatures when thesilver and iodide ions are completely uncharged, and hence no electric field ispresent (Figure 4.5).The red points and red horizontal line in Figure 4.4 are the results of Lupi et72Figure 4.4: ∆Ts required to nucleate ice in 20–200 ns versus the radius ofthe nucleating disk or plate compared with previous simulation re-sults and with the GT equation based on the properties of real water.Heterogeneous nucleation results for AgI disks or plates are shown inblack (see legend) for the six-site and TIP4P/2005 models. The solidblack line indicates the lowest supercooling required to nucleate iceon infinite AgI slabs, and the dotted black curve represents a fit ofour six-site model results to the GT equation. Red downward trian-gles and the red line show heterogeneous nucleation results[81] for themW model on graphene disks and infinite slabs, respectively. Blue tri-angles are homogeneous nucleation results for the TIP4P-E model ob-tained employing infinite cylindrical ice clusters.[66] The GT equationfor cylinders[136] is also shown (dashed blue curve). Homogeneousnucleation results employing spherical ice clusters are in green (seelegend) for the TIP4P/2005,[2] TIP4P/Ice,[2] and TIP4P[64] models.The real-water GT equation[136] for spherical nuclei is shown as agreen dashed curve.73Figure 4.5: Ice cluster size vs time for simulations done with neutral Agand I atoms (left panel) and with pressure coupling maintaining thepressure of the simulation at 1 bar. Keys show the system and ∆Ts foreach simulation.al.[81] for ice nucleation on graphene obtained using the mW water model.[6] Inthat system, ice nucleation also begins with ice-like clusters forming on the disks,and the trend observed is similar to our results; however, we note that for the mWmodel the infinite surface limit is achieved at a radius of ∼ 2 nm, whereas in ourcase the approach to the infinite limit appears to be asymptotic.In Figure 4.4 we also compare our results with predictions for homogeneousnucleation of an infinitely long cylindrical nucleus,[65, 66] homogeneous nucle-ation of a spherical critical nucleus,[2, 64] and predictions based on the Gibbs–Thomson (GT) equation, which relates ∆Ts to the size of a critical ice cluster ornucleus. The GT equation imposes an equality of the chemical potential acrossa curved interface. For small particles, which have high curvature, the surface isless stable, water molecules that are part of an ice embryo can leave easier, and74the melting point of the particle is lower than the bulk. Therefore, small clusters,having highly curved surfaces, melt significantly below the bulk melting point.The GT equation quantifies the melting point depression of a small particle vs itsbulk melting point. A second-order expansion of the GT equation can be writtenin the form[136]∆Ts =(Tmρq)γg+(κTmρq− 12χ− ∆cTm2ρq2)Tmρq(γg)2 , (4.2)where (assuming homogeneous nucleation) q is the heat of melting valued at 333J/g, ρ is the density of ice, γ is the interfacial tension between ice and water with avalue of 29 µJ/cm−2, χ is the compressibility of ice, κ is the thermal expansivity,and ∆c is the difference between the isobaric specific heats of water and ice.In equation 4.2, g depends on the geometry of the ice cluster and is 2/r for asphere of radius r and 1/r for a cylinder, where r is the radius of the base. In bothcases, equation 4.2 reduces to the form∆Ts =C1r−C2r2, (4.3)where C1 and C2 are constants that depend on the geometry. Not all of the pa-rameters determining these constants are known for the water models consideredhere, but the experimental values for water yield[136] C1 = 519 K A˚, C2 = 829K A˚2 for spherical clusters and C1 = 259 K A˚, C2 = 207 K A˚2 for the cylindricalcase.Figure 4.4 suggests that the radius of the disk or plate required for heteroge-75neous nucleation on AgI at a given supercooling agrees well with the critical ra-dius for homogeneous nucleation of spheres at the same supercooling. A similartrend has recently been observed for heterogeneous nucleation on water-solublemacromolecules.[137] In contrast, the radius of the disk or plate required for het-erogeneous nucleation on AgI is larger than the critical radius for homogeneousnucleation of infinite cylinders at the same supercooling.Three GT curves are given in Figure 4.4. The curves for spherical and cylin-drical clusters were obtained using the experimental parameters for water. Wenote that the results of Pereyra et al.[66] lie close to the GT equation for cylindri-cal nuclei. In Figure 4.4, we also include the curve for the parameters C1 = 450K A˚, and C2 = 537 K A˚2, which represent the best least-squares fit of Figure 4.3to our simulation results. We see that Figure 4.3 does give a very good fit to ourresults. Thus, the two-parameter GT equation strictly applicable to homogeneousnucleation gives a reasonably good representation of heterogeneous nucleation byAgI disks or plates, likely indicating that the ice-liquid interfacial tension remainsa dominant factor in determining critical cluster size.Figure 4.6 is a set of plots showing the number of water molecules identi-fied as ice-like forming clusters on various AgI disks and plates as time pro-gresses. Additional plots for the C2 case and for simulations with the TIP4P/2005model are shown in Figures 4.7 and 4.8, respectively. As previously described,the CHILL[130] algorithm was used to identify which water molecules are ice-like, and the DBScan[131] algorithm was employed to determine which ice-likemolecules are part of the cluster on the disk or plate surface.76Figure 4.6: Caption on next page.77Figure 4.6: Plots showing the number of water molecules identified as ice-like growing on each disk or plate at the minimum ∆Ts, where icenucleation was observed, and at a temperature one degree warmer. The∆Ts values are given in the legends. Our estimates of the minimumand maximum limits on the critical cluster size (see text) are shown ashorizontal lines and open circles, respectively.Figure 4.7: Ice cluster size versus time for the C2 case (six-site model) atthe ∆Ts values shown in the legend.Simulation results for two ∆Ts values are shown in Figure 4.6; one corre-sponds to the smallest ∆Ts for which ice nucleation was observed, and the othercorresponds to simulations one degree warmer, where ice did not nucleate. Wenote that in simulations where ice nucleation occurred, it did so in 100–200 nsafter supercooling. Of course, it is always possible that ice nucleation would oc-cur at higher temperatures (smaller ∆Ts) with longer simulation times, but oursimulations represent what is presently practical.For ice to nucleate and grow in our simulations, a critical ice cluster size (crit-78Figure 4.8: Ice cluster size versus time for the TIP4P/2005 model. The leftpanel is for disk D and the right panel for disk E. The horizontal linesindicate our lower limit estimates of the critical cluster size. Estimatesof the upper size limit are shown as circles. ∆Ts values are given in thelegend.ical nucleus) must be achieved. We estimate upper and lower limits to the numberof water molecules in the critical cluster. In Figure 4.6, lower limits are shown asblack horizontal lines. These estimates were obtained from simulations that didnot nucleate ice and that were one degree warmer than simulations that did. Thenumber of water molecules in the ice cluster was recorded every nanosecond (Fig-ure 4.6), and the average number of water molecules, together with the associatedstandard deviation, was calculated using the results for times >20 ns. The lowerlimit on cluster size was taken to be two standard deviations above the average79size obtained in these simulations.Upper limits on the number of water molecules in the critical cluster wereestimated for each simulation that nucleated ice and were taken to be the size atwhich the cluster begins to extend beyond a cylindrical volume defined by the ra-dius of the disk or plate, or, more precisely, when the number of ice-like moleculesoutside the cylinder exceeds 10. Our justification for this somewhat arbitrary se-lection of an upper limit is based on Figures 4.9, 4.10, and 4.11, where it can beseen that lateral cluster growth beyond a disk-defined cylinder closely coincideswith nucleation and irreversible growth of bulk ice. For comparison, the disk Dresults shown in Figure 4.9 are from the same simulation depicted in Figure 4.3.In Figure 4.9, the number of water molecules in the ice cluster growing on thedisk or plate is plotted as a function of time. The ice cluster is further resolvedby counting the number of molecules inside and outside the cylinder, with a basedefined by the radius of the circle circumscribing the disk or plate. Similar plotsfor other ice-nucleating simulations employing the six-site model are given inFigure 4.9, and results for the TIP4P/2005 model are shown in Figure 4.11.The plots shown in Figures 4.6, 4.9, 4.10, and 4.11 suggest that ice nucleationand growth on AgI disks and plates occur in two stages. (Note that the plotsobtained in NPT simulations (Figure 4.5) show very similar behaviour.) There isan initial period, varying in length depending on the disk or plate size, where theice cluster slowly grows, sometimes fluctuating in size, within the confines of thecylindrical space determined by the disk or plate radius, as previously described.This is followed by rapid growth once a critical cluster size is reached. Much of80Figure 4.9: Number of ice-like molecules as a function of time for four disksand plates. Curves are included for the total number of ice molecules(black) and for the numbers inside (red) and outside (green) the cylin-drical volume defined by the disk or plate radius. In each panel, a cir-cle denotes the point identified as the upper limit on the critical clustersize.the rapid growth occurs horizontally (see the green line in Figure 4.9) on whatwould be the prism face of pure Ih, where growth is known to be faster than on thebasal face[138] (perpendicular to the surface of the AgI disk or plate). There is anapparent correlation between when the cluster begins to extend beyond the disk-based cylinder and the onset of rapid growth. Therefore, we define (admittedlysomewhat arbitrarily) the transition from the first stage of growth to the secondas when the ice cluster has at least ten molecules beyond the edge of the disk or81Figure 4.10: Ice cluster size versus time for the six-site model. The totalcluster size is represented by a black curve, the amounts of ice insideand outside a cylindrical volume defined by the disk radius are redand green, respectively. The system label and the ∆Ts value are givenin the legend.82Figure 4.11: As in Figure 4.9, but for the TIP4P/2005 model.plate, and these critical clusters are marked with circles in Figures 4.6 and 4.9.The number of water molecules in the critical cluster are plotted as verticalbars joining the lower and upper limit estimates and are shown as functions of∆Ts in Figure 4.12. The CHILL algorithm is not efficient at recognizing ice-likemolecules that are immediately adjacent to a surface; therefore. we add additionalmolecules to the cluster numbers reported by CHILL and DBScan. From snap-shots of water on the surfaces, we estimate the corrections to be 21, 50, 75, 114,131, and 173 for disk or plate A, B, C, D, Dr, and E, respectively.In Figure 4.12A, we also included the number of water molecules in the criticalcluster predicted by the GT equation for homogeneous nucleation of a sphere83Figure 4.12: Caption on next page.using the parameters of real water. GT estimates for the number of ice moleculesin spherical clusters are obtained assuming an ice density of 917 kg/m3. We notethat our results do follow the trend predicted by the GT equation, but our estimateslie consistently below the GT curve. The somewhat smaller critical clusters weobserve possibly reflect the influence of heterogeneous nucleation by AgI.In panel A we also compare our estimates with previous simulations of the84Figure 4.12: Comparison of critical cluster results from our simulationswith previous experimental and theoretical work. (A) Comparisonwith previous simulations. Present results are shown as bars join-ing the upper and lower estimates for the critical cluster and are la-beled indicating the disk or plate employed. Results for the six-siteand TIP4P/2005 models are shown in black and pink, respectively.The green dotted curve represents the real-water GT equation.[136]Green points (see legend) are homogeneous nucleation results for theTIP4P/2005,[2] TIP4P/Ice,[2] and TIP4P[64] models. Blue trianglesare homogeneous nucleation results for the six-site model obtainedin an applied electric field.[77] (B) Comparison with experimentalestimates obtained using microemulsions[139] for homogeneous nu-cleation (purple squares, see legend) and heterogeneous nucleation(orange downward triangles). The open purple squares are homoge-neous nucleation results obtained employing levitated droplets.[140]Other symbols are as in panel A. (C) Simulation (present work) andexperimental results plotted versus the scaled temperature τ . Thesymbols and curves are as in panels A and B.critical cluster size for homogeneous nucleation. The results of Sanz et al.[2] em-ployed the TIP4P/2005[72] and TIP4P/ice[100] models and come from an inves-tigation of homogeneous ice nucleation using spherical Ih clusters as trial nuclei.Espinosa et al.[64] performed similar studies for the TIP4P[111] and mW[6] wa-ter models, but we only show the results for TIP4P because the results for mW aresimilar to those shown for TIP4P/ice. Yan et al.[77] investigated ice nucleation byelectric fields using the six-site model of water. Electric fields raise the meltingtemperature of water and thus allow homogeneous nucleation to be observed indirect simulations employing atomistic models but do not have a large effect onthe size of the critical nucleus for a given ∆Ts.[77] The results of Yan et al.[77](six-site model) are in good accord with the present estimates, which is perhaps85expected given that the same water model is employed in both cases. The resultsof Sanz et al.[2] (TIP4P/20005 and TIP4P/ice) and Espinosa et al.[64] (TIP4P) arehigher than our estimates for heterogeneous nucleation on AgI surfaces. This dif-ference may partially be due to the different models considered or to the fact thatSanz et al.[2] and Espinosa et al.[64] used spherical ice clusters that will presenta variety of crystal faces for nucleation. Different faces will grow differently andmay impede nucleation, thus requiring larger sizes. In general, panel A suggeststhat the number of water molecules in the critical cluster for heterogeneous nucle-ation on AgI is within a factor of 4 of the number of water molecules in the criticalcluster required for homogeneous nucleation of spheres at the same supercooling.To test for possible water model sensitivity, we performed simulations withdisks D and E using the TIP4P/2005 model considered by Sanz et al.[2] We findthat the TIP4P/2005 model does require a larger ∆Ts (compared with the six-sitemodel) to nucleate on disks D and E, while giving critical clusters similar in sizeto those found for the six-site model (see Figure 4.12A,C). Thus, the particularwater models employed to study ice nucleation are also a factor in determiningthe critical cluster size.Some experimental results for the critical cluster of ice obtained under differ-ent conditions are compared with our estimates in Figure 4.12B. Liu et al.[139]investigated both homogeneous and heterogeneous nucleation using water in oilmicroemulsions. To investigate heterogeneous nucleation they incorporated hep-tacosanol as the ice nucleating agent into the emulsions. In addition, Kramer etal.[140] investigated homogeneous nucleation using levitated single water droplets.86Our results show reasonable agreement with the homogeneous studies, which isconsistent with the conclusions reached from Figure 4.12A. For the case of hetero-geneous nucleation on heptacosanol, our results deviate from the measurementsat warmer temperatures, which is discussed later.Malaspina et al.[141] have suggested that another way of comparing differentwater models is to plot them versus a scaled temperature, τ , as is done in Figure4.12C. The scaled temperature is defined as τ = (T − Tw)/(Tm− Tw), where Tis the nucleation temperature, Tm is the melting temperature of ice for a particu-lar model, and Tw is the temperature of the Widom line, defined by heat capacitymaxima. The Widom line temperatures at 1 atm were taken from Malaspina etal.[141] to be 249 K for six-site, 229 K for TIP4P/2005, and 235 K for real wa-ter. Widom line temperatures are not known for the TIP4P and TIP4P/Ice models,and these models are not included in panel C. The melting temperatures used are289 K for six-site[113] and 250.5 K for TIP4P/2005.[112] The scaled temperatureis an attempt to account for dynamical differences between real water and watermodels and also between different water models that often have significantly dif-ferent melting points[141] and dynamical properties. We note that on the τ scalethe TIP4P/2005 results (both those of Sanz et al.[2] and those we obtain) stilldiffer significantly from the estimates for the six-site model. Thus, model differ-ences that are not captured by τ clearly persist and should be kept in mind in icenucleation studies of water models.From Figure 4.12B,C, we see that our model estimates and the experimentalresults for heterogeneous nucleation on heptacosanol are in reasonable agreement87at lower temperatures but deviate at the higher temperatures. At higher temper-atures, the experiments report smaller amounts of water in critical clusters andappear to follow a different trend. The higher temperature results represent theonset of heterogeneous nucleation using the IN heptacosanol. Our results are alsofor heterogeneous nucleation but appear to more closely follow the trend of ho-mogeneous nucleation.The apparent similarity of heterogeneous nucleation by AgI disks to the ho-mogeneous case in terms of the critical cluster size may reflect the fact that AgIis a close ice mimic (Chapter 3). In other words, from the point of view of icenucleation and growth, there may be a strong resemblance between a AgI disk orplate and a corresponding piece of ice, the crucial difference being that the AgIdisk or plate obviously cannot melt.Some additional insight can be gained by examining the sizes of subcriticalclusters that never achieve nucleation. We collect samples of these failed clustersusing the local maxima in Figure 4.6 and in similar plots (not shown) at othertemperatures. An example plot illustrating how we identified failed clusters forplate C at ∆Ts = 19 K is given in Figure 4.13. Note that only the larger failedclusters are selected.The number of water molecules in failed clusters for different disks or plates(denoted with different colours) is plotted versus ∆Ts in Figure 4.14. The numbersof water molecules in critical nuclei are included for comparison. Two trends areimmediately obvious from Figure 4.14. First, as one would expect, larger disksand plates create larger clusters. Second, for a given disk or plate, the number of88Figure 4.13: An example of the peaks selected for the failed cluster plot(Figure 7). This simulation was for the six-site model with plate C at∆Ts = 19 K.water molecules in the largest failed cluster decreases with increasing temperature.This second observation indicates an additional hurdle to ice nucleation at warmertemperatures. Not only are larger clusters required to nucleate ice at warmer tem-peratures but also it becomes harder to create larger clusters as the temperature isincreased. Therefore, in our simulations the warmest temperature where a partic-ular AgI disk or plate can nucleate bulk ice is the temperature where the requirednumber of water molecules in the critical cluster just matches the largest clusterthat can form on the disk or plate.89Figure 4.14: Sizes of subcritical clusters (six-site model) that failed to nucle-ate for different disks and plates as functions of ∆Ts. Critical clustersize ranges for the smallest ∆Ts where nucleation occurred are shownfor each disk or plate as bars joining estimates of upper and lowerlimits.4.6 ConclusionsWe have used molecular dynamics simulations to investigate ice nucleation by AgIdisks or plates with circumscribing radii ranging from 1.15 to 2.99 nm. The six-site and TIP4P/2005 water models were considered, and a similar ice nucleationprocess was observed in both cases. The disk or plate size required to nucleate iceon simulation time scales increases as the degree of supercooling ∆Ts decreases.For the six-site model the largest disk considered (radius of 2.99 nm) nucleatedice at ∆Ts = 14 K, compared with 8 K for an infinite AgI sheet (Chapter 3). Ourresults demonstrate that relatively small patches of defect-free AgI surface couldserve as efficient active sites for ice nucleation.Ice nucleation by nanoscale AgI disks and plates occurs in two stages. Ini-90tially, upon supercooling, ice clusters roughly hemispherical in shape form on thesilver exposed face of the AgI disk or plate. These clusters fluctuate in size butnever melt away until (if the disk size and ∆Ts are suitably matched) the clusterachieves critical size, allowing bulk ice to nucleate and rapidly grow to fill the sim-ulation cell. For a given ∆Ts, our estimates of the critical cluster size are roughlycomparable to those obtained from previous simulation studies of homogeneousnucleation. Additionally, the critical cluster size as a function of ∆Ts is well repre-sented by a second-order (two-parameter) GibbsThomson equation, as one wouldexpect for homogeneous nucleation. Thus, it appears that the main influence ofthe AgI disk or plate is to create and stabilize an ice cluster of some minimum size(dependent on ∆Ts and the radius of the AgI disks or plate) that is always presentin the system. This constantly present ice cluster reduces the free-energy barrierthat must be overcome to achieve critical cluster size and hence ice nucleation.91Chapter 5Simulations of Ice Nucleation byKaolinite (001) with Rigid andFlexible Surfaces“Beware of the man who works hard to learn something, learns it,and finds himself no wiser than before,” Bokonon tells us. “He is fullof murderous resentment of people who are ignorant without havingcome by their ignorance the hard way.”— Kurt Vonnegut, Cat’sCradle5.1 SummaryNucleation of ice by airborne particles is a process vital to weather and climate, yetour understanding of the mechanisms underlying this process is limited. Kaoliniteis a clay that is a significant component of airborne particles and is an effective icenucleus. Despite receiving considerable attention, the microscopic mechanism(s)92by which kaolinite nucleates ice is not known. We report molecular dynamicssimulations of heterogeneous ice nucleation by kaolinite (001) surfaces. Both theAl-surface and the Si-surface nucleate ice. For the Al-surface, reorientation ofthe surface hydroxyl groups is essential for ice nucleation. This flexibility allowsthe Al-surface to adopt a structure which is compatible with hexagonal ice, Ih, atthe atomic level. On the rigid Si-surface, ice nucleates via an unusual structurethat consists of an ordered arrangement of hexagonal and cubic ice layers, joinedat their basal planes where the interfacial energy cost is low. This ice structureprovides a good match to the atomistic structure of the Si-surface. This exampleis important and may have far-reaching implications because it demonstrates thatpotential ice nuclei need not be good atomic-level matches to particular planes ofice Ih or cubic ice, Ic. It suggests that surfaces can act as effective ice nuclei bymatching one of the much larger set of planes that can be constructed by regulararrangements of hexagonal and cubic ice.5.2 MethodsFor kaolinite simulations in Chapter 5, several variations of the Al-surface areconsidered: (1) Simulations with kaolinite held completely fixed with all atoms intheir crystallographic positions are referred to as rigid-kaolinite. (2) The hydro-gen atoms of the Al-surface hydroxyl groups are allowed to move, but they remainbonded to their oxygen atoms with a fixed bond length. This is referred to as thefree-H surface. These conditions give the effect of the hydroxyl group pivotingon its oxygen atom such that the OH bond adopts various orientations. (3) In93addition to the free-H conditions, we also sometimes allow the hydroxyl oxygenatoms to vibrate within a harmonic potential of 1000 kJmol−1nm−2, similar to thepotential employed by Fraux and Doye[83] for silver iodide. This is referred to asthe free-OH surface. We attempted to simulate a layer of kaolinite with all atomscompletely free to move by supporting it from below with two fixed layers,[142]but the free layer failed to remain intact. Allowing the hydroxyl oxygen atomsto vibrate imparts some of the motion expected in an entirely free surface, with-out the risk of the surface falling apart. In the free-OH case, the OH bond lengthis again held fixed. (4) Lastly, several rigid Al-surface structures with noncrys-tallographic, or prepared, hydrogen atom coordinates are considered. We alsosimulate the Si-surface using crystallographic coordinates. The Si-surface doesnot have any hydroxyl groups on it, and in our simulations we do not add anyhydroxyl groups. Dimensions of the simulations used for kaolinite are given inTable 5.1. Lists of simulations used for Chapter 5 are provided in Tables C.1–C.9.Table 5.1: For each system considered, the XY cell dimensions, the numberof water molecules, and the bulk water densities at 300 K are given.The Z cell dimension is 8.42661 nm in all cases (surfaces are 6 nmapart measured from the edges of the unit cells), except for the largeSi-surface, which is 1 nm larger.System X (nm) Y (nm) Water molecules Density (g/mL)Al-surface 5.66885 6.2593 7100 0.940-0.943Si-surface 5.66885 6.2593 6800 0.944large Si-surface 10.307 10.73016 17700 0.930945.3 Results5.3.1 Varying Flexibility of the Al-SurfaceFigure 5.1 shows ice growth as a function of time in several simulations for Al-surfaces with varying degrees of flexibility. Initially, there is little ice in the free-Hsimulation, but ice begins to form after 75–250 ns and continues to grow until thesimulation cell is full, and the plots level off. The final state of a free-H simulationwith six-site water is shown in Figure 5.2. Ice Ih is the dominant form of iceformed in these simulations, and the figure is oriented so that one is looking alongthe c-axis (into the basal plane) of ice Ih. A few layers of ice Ic also formed inthis simulation but do not make contact with the surface, and they are the watermolecules visible within the holes of the ice Ih lattice.95Figure 5.1: Number of ice molecules as a function of time for three differ-ent Al-surfaces, as indicated in the legends. Results are shown forthe TIP4P/Ice (top) and six-site (bottom) water models. The rigid-kaolinite, free-H, and free-OH surfaces are described in the text.The other simulations shown in Figure 5.1 did not result in ice nucleation.These include the rigid-kaolinite case, in which the kaolinite hydrogen atoms arefixed onto their crystallographic coordinates and are not able to adapt to the pres-ence of water molecules and the free-OH surface. In the free-OH simulations,the hydroxyl hydrogen atoms are free to move; as in the free-H case, the oxygen96Figure 5.2: Configurational snapshot of ice grown on the Al-surface (free-H) with the six-site model. The simulation cell is oriented such thatone is looking along the c-axis (into the basal plane) of ice Ih. Watermolecules visible within the hexagons formed by the ice Ih lattice arefrom a few sheets of ice Ic. Mirrored kaolinite slabs are at the topand bottom of the simulation cell. Hydrogen, aluminum, silicon, andoxygen atoms are black, pink, yellow, and red, respectively, except foroxygen atoms reported as ice, which are blue.atoms are also free to move, but are restrained to their crystallographic positionsby means of a harmonic potential. The free-OH construction is meant to ap-proximate an unconstrained kaolinite surface while avoiding the possibility of thesurface layer breaking apart (atoms leaving the surface, etc.) during the simula-tion. However, as noted above (also see discussion below), we did not observe icenucleation on a free-OH surface.97The three simulations (rigid-kaolinite, free-H, and free-OH) reveal that flex-ibility (but not too much) of the kaolinite surface results in ice nucleation. Thestructure of water that forms on the crystallographic Al-surface is not compatiblewith ice; hence, nucleation does not occur in the rigid-kaolinite case, consistentwith earlier simulation results.[87, 88] On the other hand, if the surface atoms areallowed too much freedom, as in the free-OH case, we did not observe ice nucle-ation on simulation time scales. This might be what prevented ice nucleation inearlier work.[86]We carried out a set of simulations to investigate if the free-OH surface issomehow incompatible with ice, or if some kinetic aspects of nucleation are in-hibited by motion of the oxygen atoms. These simulations were initiated in free-Hmode; then, at a later time they were switched to free-OH conditions (hydroxyloxygen atoms are allowed to vibrate about their lattice positions as describedabove) and continued to see if ice formed during the initial free-H period wouldcontinue to grow, or possibly melt. The results (number of ice molecules versustime) are shown in Figure 5.3. The free-H portions of the simulations are shown ingray, and the free-OH parts as coloured curves. Configurational snapshots show-ing ice growth (or not) for these simulations are shown in Figures 5.4 and 5.5.We note that if some ice forms (approximately 2 surface layers or 250–500 icemolecules are sufficient) during the free-H period of the simulation, then ice con-tinues to grow after switching to free-OH conditions. This suggests that ice iscompatible with the free-OH surface, but that the vibration of the hydroxyl oxy-gen atoms increases the time required for nucleation beyond the time scale of our98Figure 5.3: Number of ice molecules as a function of time on Al-surfacesfor simulations begun under free-H conditions and then switched tofree-OH. The parent free-H parts are shown as gray curves, and thefree-OH continuations are indicated in colour. Results are shown forthe TIP4P/Ice (top) and six-site (bottom) water models.simulations.Further insight into the water (or ice) structure near the three Al-surface varia-tions can be obtained from the density profiles and snapshots presented in Figure5.6 for the TIP4P/Ice model. Similar results for the six-site water model are given99Figure 5.4: Caption on next page.100Figure 5.4: Snapshots showing the beginning and end of each of the free-OHsimulations begun with configurations taken from free-H simulationsfor the six-site model. Each pair of snapshots shows the beginning (topimage) and end (bottom image) of the free-OH simulation. The firsttime given above each pair of snapshots is the length of the initial free-H simulation, and the second time is the length of the simulation underfree-OH conditions. The colours for oxygen, hydrogen, aluminum,and silicon atoms are red, black, purple, and yellow, respectively.Figure 5.5: Snapshots showing the beginning and end of each of the free-OHsimulations begun with configurations taken from free-H simulationsfor the TIP4P/Ice model. Format and colours are as in Figure 5.4.101in Figure 5.7. The four panels in Figure 5.6 show results for rigid-kaolinite (5.6A),free-OH (5.6B), and free-H before (5.6C) and after (5.6D) ice nucleation. Den-sity profiles for the oxygen and hydrogen atoms of water, and in some plots thoseof the kaolinite hydroxyl groups, are shown. The zero point of the y-axis of thedensity plots is set at the crystallographic position of the kaolinite oxygen atoms.A configurational snapshot is positioned adjacent to each density profile, and asnapshot of the first water layer on the particular kaolinite surface is shown beloweach profile.First we consider Figure 5.6D which shows the density profiles for ice thathas grown on the free-H surface. We note that the profile for the water oxygenatoms is a series of doublets, reflecting the bilayers that comprise the hexagonalrings of ice. The water hydrogen atom profile has strong peaks lying betweenthe oxygen atom doublets that are due to hydrogen atoms hydrogen bonded to theoxygen atoms. On both sides of the oxygen doublets are weaker hydrogen peaksformed by hydrogen atoms that hydrogen bond to the kaolinite surface and/or tobilayers above and below. The density profile of the kaolinite hydroxyl hydrogenatoms shows that they are usually directed away from the surface, but ∼ 1/4 lieparallel to the surface. The snapshot below the density profile shows that the watermolecules form hexagonal rings on the surface. Results for liquid water near thefree-H surface are shown in panel C. We note that the surface imprints icelikestructure in the nearby liquid. This is similar to our earlier observations for somesurfaces of silver iodide, which also induced icelike structure in the liquid nearthe surface, prior to bulk freezing (Chapter 3).102Figure 5.6: Four panels illustrating the behaviour of TIP4P/Ice water as wellas the kaolinite hydroxyl groups near Al-surfaces. The surfaces in-cluded are rigid-kaolinite (A), free-OH (B), free-H before ice nucle-ation (C), and free-H after ice nucleation and growth (D). Each panelcontains density profiles (top left), a configurational snapshot showingwater and surface hydroxyl groups (top right), and a snapshot of thefirst layer of water near the Al-surface and the surface hydroxyl groups(bottom). Particular atoms are indicated with the same colour in alldensity profiles and snapshots: water hydrogen atoms (black), wateroxygen atoms (red), hydroxyl hydrogen atoms (blue), and hydroxyloxygen atoms (green). Density profiles are aligned with the snapshotsto their right, and the zero on the y-axis coincides with the crystal-lographic position of the hydroxyl oxygen atoms of the Al-surface.Density profiles are shown only for atoms that can move during thesimulation. Therefore, panel A includes only density profiles for theoxygen and hydrogen atoms of water, and the density profile of thekaolinite hydroxyl oxygen atoms is shown only in panel B. The yellowhexagons in panels C and D (bottom) highlight examples of hexagonalrings formed by water in the surface layer.103Figure 5.7: Four panels illustrating the behaviour of six-site water as well asthe kaolinite hydroxyl groups near Al-surfaces. The surfaces includedare rigid-kaolinite (A), free-OH (B), free-H before ice nucleation (C),and free-H after ice nucleation and growth (D). Each panel containsdensity profiles (top left), a configurational snapshot showing waterand surface hydroxyl groups (top right), and a snapshot of the firstlayer of water near the Al-surface and the surface hydroxyl groups(bottom). Particular atoms are indicated with the same colour in alldensity profiles and snapshots: water hydrogen atoms (black), wateroxygen atoms (red), hydroxyl hydrogen atoms (blue), and hydroxyloxygen atoms (green). Density profiles are aligned with the snapshotsto their right, and the zero on the y-axis coincides with the crystal-lographic position of the hydroxyl oxygen atoms of the Al-surface.Density profiles are shown only for atoms that can move during thesimulation. Therefore, Panel A only includes density profiles for theoxygen and hydrogen atoms of water, and the density profile of thekaolinite hydroxyl oxygen atoms is shown only in Panel B.104The free-OH profile (panel B) is very similar to the free-H liquid profile shownin panel C and also resembles results of a previous density functional theory (DFT)study,[51] suggesting that the force fields employed here give a reasonable pictureof the behaviour of the kaolinite surface. The strong similarity of the profiles inpanels B and C suggests that the free-OH surface (panel B) is compatible with iceand, as discussed above, longer simulations might lead to ice nucleation. Notethat the kaolinite hydroxyl oxygen atoms (green curve in the density profile) haveonly a narrow range of oscillation.The rigid-kaolinite profile (panel A) has some similarity to panel C (free-H),but there are apparently crucial differences. In panel C, the first water oxygenpeak lines up with the leading edge of the second water hydrogen peak, whereasthis is not the case for rigid kaolinite. Additionally, for rigid-kaolinite (panel A)the first hydrogen peak lies closer to the surface. These differences may appearsmall at first sight, but they indicate a less icelike structure in the adjacent water.This is likely why ice nucleation was not observed for the rigid-kaolinite surfacein the present and earlier simulations.[87–90]Figure 5.8 shows the amounts of ice Ih and ice Ic formed during simulationsemploying the free-H surface. We see that for both the six-site and TIP4P/Icemodels, hexagonal ice is by far the dominant structure. For the six-site model oneor two layers of cubic ice formed toward the end of the simulation, but we note(Figure 5.2) that these cubic layers do not make contact with the kaolinite surface.Calculation of the lattice mismatch clearly shows that the Al-surface of kaoli-nite is expected to favor ice Ih. The lattice mismatch compares an ice lattice with105Figure 5.8: Number of ice Ih and ice Ic molecules as a function of time grow-ing from the Al-surface under free-H conditions. Results are shown forthe six-site (240 K) and TIP4P/Ice (230 K) water models, as indicatedin the legend. Note that ice Ih is by far the dominant ice form.that of a potential heterogeneous ice nucleus via the equationmismatch =Inan− ImaiImai100% , (5.1)where an and ai are the lattice constants of the ice nucleus and of ice, respectively,and In and Im are integers chosen to minimize the mismatch.[1] Mismatch per-centages are given in Table 5.2. A good match between kaolinite and ice will haveice matching kaolinite in both directions of the unit cell, with mismatches of onlya few percent. The results given in the table show that kaolinite (001) is a goodmatch only for the prism face of ice Ih (5.7%, 0.2%), consistent with what we106Table 5.2: Lattice mismatch percentages from equation 5.1. Numbers insquare brackets give the values of Im and In ([Im : In]) used. Numbersin parentheses are the lattice parameters. The only good match is high-lighted in bold.lattice parameter Kaolinite-a (5.15 A˚)[118] Kaolinite-b (8.94 A˚)[118]Ih-a (4.48 A˚)[18] [1:1]15 % [2:1]0.2 %Ih-c (7.31 A˚)[18] [2:3]5.7 % [1:1]22.3 %Ic (6.37 A˚)[18] [1:1]19.1 % [3:2]6.4 %see in the simulations. There is no good match for ice Ic; hence, we see little orno cubic structure in ice nucleated on the free-H surface. This result agrees withprevious simulation work of Cox et al.[86] who observed that ice nucleating incontact with a kaolinite surface was exclusively ice Ih. The results given in Table5.2 were obtained using experimental values for the unit cells of ice. Ice Ih unitcell values for six-site water have been reported[110] and give mismatches within12% of those in Table 5.2.It is important to note that while a good lattice match explains why the prismface of ice Ih nucleates on kaolinite (001), it is the reorientation of the surfacehydroxyl groups that allows ice to nucleate at all. The hydrophilic[50] hydroxidesurface of kaolinite (001) is covered with hydroxyl groups that take one of twoprimary orientations. Hydroxyl groups can point away from the surface and do-nate hydrogen bonds to water, or lie parallel to the surface, such that the oxygenatoms can accept a hydrogen bond from water. This behaviour allows the Al-surface to be amphoteric and led previous workers[93–95] to ascribe kaolinitesice-nucleating ability to this behaviour.1075.3.2 Non-crystallographic Rigid Al-surfacesFurther insight into the role of the orientation of the surface hydroxyl groups inice nucleation can be gained by examining several rigid surface structures pre-pared in different ways. We consider four such rigid surfaces. Two were obtainedby taking surface configurations from free-H simulations (at 230 K for TIP4P/Iceand 240 K for six-site) before and after freezing, which we label liquid-preparedand ice-prepared, respectively. A third ideal surface was constructed based on ourobservations for the ice-prepared surface, and a fourth surface was obtained byvacuum annealing (no water present) the kaolinite surface at 25 K. The surfaceproduced by vacuum annealing closely resembles the surfaces predicted by DFTsimulations.[93–95] Examples of these prepared rigid surfaces are shown in Fig-ures 5.9, 5.10, and 5.11. In these figures, only the oxygen atoms of the surfacehydroxyl groups and the hydroxyl hydrogen atoms that are directed away from thesurface are shown. Also, the hydrogen atoms are enlarged such that they obscurethe underlying oxygen atoms. This surface view allows one to easily observe anypattern formed by the orientations of the hydroxyl groups.The results obtained for the four surfaces are summarized in Table 5.3. Ofthe four surfaces considered, only the ice-prepared surface and the ideal surfacenucleated ice in our simulations. The ice-prepared surface nucleated ice immedi-ately upon cooling and did so up to 255 K (15 K below the melting temperature)for the TIP4P/Ice model and up to 260 K (29 K below the melting temperature)for six-site water. Ice nucleation at these temperatures also occurred when ice-prepared configurations from one water model were used with the other model,108Figure 5.9: Hydroxyl group configurations used in the prepared rigid surfacesimulations. Oxygen atoms are red, and hydrogen atoms are blue; onlyhydrogen atoms directed upward are shown.showing that the favorable surface configurations are not model-specific. Thefact that the rigid ice-prepared surface readily nucleates ice suggests that it is thesurface structure, rather than some dynamical effect related to the motion of sur-face hydrogen atoms, that is important for ice nucleation. Once an appropriatestructure is achieved, nucleation occurs whether or not the hydroxyl groups areallowed to reorient during the simulation. The liquid-prepared surface bears someresemblance to the ice-prepared case (Figure 5.9) but did not nucleate ice in oursimulations. This suggests that the final surface structure complementary to ice isestablished during the freezing process.109Figure 5.10: Additional ice- and liquid-prepared surfaces obtained for theTIP4P/Ice model. Hydroxyl oxygen atoms are red, hydroxyl hydro-gen atoms are blue, and only hydrogen atoms directed upward areshown.110Figure 5.11: Additional ice- and liquid-prepared surfaces obtained for thesix-site model. Hydroxyl oxygen atoms are red, hydroxyl hydrogenatoms are blue, and only hydrogen atoms directed upward are shown.111Table 5.3: Prepared rigid surfaces. Yes and No indicate if ice nucleated ornot.Preparation six-site TIP4P/Iceideal Yes Yesvacuum No Noice (same model) Yes Yesice (different model) Yes Yesliquid No NoUpon close examination, ice-prepared surfaces appear to have a little moreorder than liquid-prepared surfaces. In particular, ice-prepared surfaces have arepeating pattern of clusters of upward pointing hydrogen atoms, highlighted withyellow hexagons in Figure 5.9. These clusters consist of a hexagonal ring ofsix hydrogen atoms surrounding a central hydrogen atom. Despite some obviousdefects, these clusters appear to be regularly positioned to form a larger hexagonallattice. On the basis of these observations, we constructed the ideal ice nucleationsurface shown in Figure 5.9. This surface nucleated ice Ih, but nucleation requiredmore time (∼ 200 ns) than ice-prepared surfaces. Thus, the pattern of the idealsurface clearly traps the essential features required for ice nucleation on the Al-surface, but the slower nucleation time suggests that some finer details might alsofacilitate ice nucleation.The final rigid surface considered is the vacuum-prepared surface shown inFigure 5.9. As noted above, this surface closely resembles the structure ob-tained in DFT calculations[93–95] (providing some confirmation of the accuracyof the CLAYFF force field[122]), but it did not nucleate ice. We note that for thevacuum-prepared surface 1/3 of the hydroxyl hydrogen atoms are lying parallel to112the surface, compared to ∼ 1/4 in the ice-prepared case.5.3.3 Si-surfaceAs previously mentioned, kaolinite (001) has two different faces, both with thesame crystallographic match to ice, but with different composition and structureat the atomic level. We also carried out simulations attempting to nucleate iceon the hydrophobic[50] Si-surface shown in Figure 5.12A. Unlike the Al-surfacediscussed above, the Si-surface was held rigid in all simulations. All simulationscarried out with the Si-surface are summarized in Table C.9. For the Si-surface,we observed ice nucleation with both water models. Ice nucleation on the Si-surface is interesting in that it followed a path different from any that we hadpreviously observed. Specifically, the first layers of ice to form on this surface donot correspond to any bulk phase of ice of which we are aware. This result wasrepeated in several simulations with both models.A series of snapshots showing ice growth on the Si-surface is shown in Figure5.13. Initially, an ordered layer of water forms on the surface composed of five-, six-, and eight-membered rings. Examples of six- and eight-membered ringsare highlighted in Figure 5.12B; five-membered rings stand perpendicular to thesurface and can be seen in Figures 5.13 and 5.14.Figure 5.13 shows that unlike the Al-surface, the basal planes of ice are notstacked perpendicular to the Si-surface. Also, unlike the single hexagonal icephase that nucleates on the Al-surface, for the Si-surface the CHILL algorithmdetects layers of both hexagonal and cubic ice. However, this is not a random113ABFigure 5.12: Caption on next page.114Figure 5.12: Si-surface without (A) and with (B) a layer of water. The oxy-gen atoms of water are blue, and hydrogen atoms are black. Theoxygen atoms of kaolinite are red, and silicon atoms are yellow. Theyellow and green lines in panel B outline the eight- and six-memberedrings, respectively, formed by water on the Si-surface.Figure 5.13: Configurational snapshots of ice nucleating and growing on theSi-surface. The results shown are for the six-site model at 240 K. Thecolours are as in Figure 5.2 except that oxygen atoms of ice Ih andice Ic are blue and green, respectively.115(0110)(0331)Figure 5.14: Close-up of the snapshot at 700 ns from Figure 5.13. The redarrows and indices indicate planes in ice Ih, and the blue arrow tracesthe layer of ice in contact with the Si-surface.mixture of hexagonal and cubic ice layers, such as was observed in previous ex-periments and simulations[24–26] and known as ice Isd . Rather, in the presentcase, the ice that grows is adapted to fit the surface structure with ordered layersof hexagonal and cubic ice (as reported by CHILL). Close inspection of snapshots(Figure 5.14), with the aid of a molecular model (Figure 5.15), reveals that the116Figure 5.15: A photograph of a model of the ice formed on the Si-surface.Black balls represent water oxygen atoms. The Si-surface is repre-sented by the black background at the bottom of the picture.stacking of hexagonal (h) and cubic (c) layers follows the pattern hhhhcchhhhcc· · · .Also, we note that the basal planes lie at an angle of 13−15◦, with respect to thesurface normal.A view of the structure is given in Figure 5.14. This is a close-up of the 700ns configurational snapshot shown in Figure 5.13, highlighting the orientationof the hexagonal planes. If we consider the ice Ih structure, prism planes areperpendicular to the basal plane, and are denoted by (011¯0) in Figure 5.14. The(033¯1¯) plane lies at an angle to (011¯0) and is parallel to the Si-surface. If wefollow the red arrow indicating the (033¯1¯) plane, we note that it passes throughthe bottom of two hexagonal rings (as we would expect for ice Ih), but throughthe middle of the next ring, which lies at the junction of two sets of basal planes,and is a ring from ice Ic. The blue arrow parallel to the surface passes though the117bottom of two hexagonal rings, followed by the bottom of a five-membered ring.The five-membered ring is a hexagonal ring disrupted by the surface. Thus, thelayer of ice in contact with the Si-surface bears some resemblance to the (033¯1¯)plane ice Ih but contains regular defects.To be sure that our Si-surface results are not strongly influenced by finite-sizeand/or the fixed cell volume, we carried out a simulation of one system approxi-mately four times larger and with a free surface essentially corresponding to zeropressure (see large Si-surface in Table 5.1 and Figure 5.16). The system was equi-librated at 300 K for 1 ns, after which the top slab was pulled back 1 nm fromthe water (with a corresponding 1 nm extension of the simulation cell), and thesystem was cooled to 250 K. Ice nucleation proceeded in a manner similar to theother Si-surface simulations, verifying the absence of any significant dependenceon system size or confinement.As noted above, earlier work has shown that mixtures of hexagonal and cubicice occur quite commonly, which is not surprising as both phases of ice containbasal planes that bond to each other with only a small interfacial energy.[23] Thedifference here is that on the Si-surface the stacking pattern of hexagonal andcubic layers is controlled by the surface to produce a regular structure. As a bulkphase, we believe that the ice formed on the Si-surface would have a monoclinicunit cell.An interesting outcome of the Si-surface result is the realization that the num-ber of ice planes that can match a particular surface is not limited to the planesof ice Ih or ice Ic, but extends to what must be a much larger number of potential118Figure 5.16: Snapshots of ice nucleating (six-site model at 250 K) on thelarge Si-surface. Hydrogen, aluminum, silicon, and oxygen atoms areblack, pink, yellow, and red, respectively, except for oxygen atoms ofice, which are blue for ice Ih and green for ice Ic.matches, given the myriad of possible stacking patterns involving combinations ofhexagonal and cubic layers. This appears to be a hitherto unrecognized possibilityfor heterogeneous ice nucleation.5.4 ConclusionsIn this paper we report simulations of heterogeneous ice nucleation on two kaoli-nite (001) surfaces. The simulations reveal interestingly different ice nucleationprocesses for each face. Surface adaptation of the Al-surface allows nucleation ofice Ih, whereas nucleation occurs via a unique surface-adapted ice structure on the119Si-surface.Ice nucleation on the Al-surface did not occur with the hydroxyl hydrogenatoms fixed in their crystallographic positions. Ice nucleation on this surface re-quired that the surface hydroxyl groups reorient to adopt a configuration compati-ble with ice. In the configuration favorable to ice nucleation, the surface hydroxylgroups act as both donors and acceptors in hydrogen bonds with water moleculesin the ice lattice. On the Al-surface, ice Ih is the dominant ice structure observed,consistent with expectations based on lattice mismatch calculations. We speculatethat the adaptive surface mechanism may be relevant to ice nucleation by otherflexible ice nuclei.We also obtained ice nucleation on the rigid Si-surface. Ice nucleation at thissurface occurs by a different and interesting mechanism. The atomic structure ofthe Si-surface is not a good atomistic match to any plane of ice Ih or ice Ic, and apure ice Ih or ice Ic phase does not grow at this surface. Rather, one observes anice structure in which layers of hexagonal and cubic ice form a repeating patternat a particular angle to the Si-surface. This ice structure is imposed by the detailedatomic structure of the Si-surface. The ice formed is similar to ice Isd in that thehexagonal and cubic layers join at basal planes where the interfacial energy costis low. However, it is unlike ice Isd in that the stacking is not disordered but formsa regular pattern dictated by the surface structure. This mechanism of nucleationopens the possibility for ice nucleation on surfaces that may be good atomic-levelmatches to a variety of planes which do not exist in either ice Ih or ice Ic but canbe found in some regular combination of these two structures.120Chapter 6Ice Nucleation Abilities of PotassiumFeldspar and Gibbsite“No wonder kids grow up crazy. A cat’s cradle is nothing but a bunchof X’s between somebody’s hands, and little kids look and look andlook at all those X’s...”“And?”“No damn cat, and no damn cradle.”— Kurt Vonnegut, Cat’s Cradle6.1 SummaryAtmospheric aerosols, which nucleate ice, are an important component of earth’sweather, initiating a large portion of rain. Clays and minerals are a large source ofatmospheric aerosols and ice nuclei. Potassium feldspar has recently been discov-ered to be a very effective ice nucleus in experimental freezing studies.[46, 52]We performed molecular dynamics simulations on the (001) and (010) faces ofpotassium feldspar and did not observe any ice nucleation. We recognize that121experimental samples, and aerosol particles, of clays and minerals can containmultiple species, so we also investigated gibbsite, which is a weathering productof feldspar. Our simulations were not able to directly nucleate ice on gibbsite.However, we were able to produce a preconfigured gibbsite surface which was aneffective ice nucleus. From these results we speculate that weathering products offeldspar may account for some of feldspar’s ice nucleating ability.6.2 MethodsFeldspar contains two primary cleavage planes,[49] the (001) and (010) planes.We cleave these planes such that oxygen atoms lying on these planes all remainattached to one of the surfaces. Keeping all the oxygen atoms attached to oneside of the cut produces two faces: a bare surface whose outermost atoms aresilicon, aluminium, and potassium, and the +O surface, which is the bare surfacewith the bridging oxygen atoms attached. A third surface is then constructed byprotonating the oxygen atoms, which we refer to as the +OH surface. See Figure6.1 for images of each of the surfaces. The hydroxyl groups were treated as in thekaolinite free-H simulations, the hydrogen atoms were allowed to pivot around thefixed oxygen atom to which they were attached. Feldspar cleaves smoothly, andhydroxyl groups form on the surface,[143, 144] so our simulations should cover arange of realistic possibilities.The CLAYFF force field[122] provides interaction parameters for feldspar.For the bare and +O surfaces, all the partial charges balance, however, addition ofhydrogen atoms for the +OH surface creates a net charge which must be neutral-122+OH     +O      bare(010)    (001)Figure 6.1: Images of potassium feldspar surfaces considered. Silicon, alu-minium, potassium and bridging oxygen atoms are yellow, teal, pink,and red respectively. Hydroxyl groups, and surface oxygen atoms inthe +OH and +O surfaces are purple. The top row displays the (001)surface and the bottom row displays the (010) surface. The type ofsurface is indicated at the top of each column.ized. Our solution for the +OH surfaces, was to add small amounts of negativecharge to all the atoms in feldspar, except potassium, to balance the extra charge.Different amounts of charge were given to each type of atom, but in all cases, thecharges varied by less than 4.5%. The net charge of the (010)+OH simulationcould not be fully neutralized by changing partial charges. To balance the chargesin this simulation, a single point charge of -0.0011e was placed in the vacuumportion of the simulation cell equidistant from each feldspar slab. Simulations onfeldspar were run at temperatures of either 240 or 250 K for six-site water, andat 230 or 240 K for TIP4P/Ice. Summaries of the simulation cells explored are123provided in Table 6.1.Table 6.1: Simulation parameters and system properties. For each systemconsidered, the XY cell dimensions, the number of water molecules,and the bulk water densities at 300 K are given.Surface X (nm) Y (nm) Water molecules Density (g/mL)Feldspar(001)+OH 5.99424 6.4815 7360 0.951Feldspar(001)+O 5.13792 5.1852 4400 0.984Feldspar(001) 5.1379 5.1852 4400 0.943Feldspar(010)+OH 5.1378 5.1808 4140 0.932Feldspar(010)+O 5.1378 5.1808 4400 1.007Feldspar(010) 5.1378 5.1808 4400 0.941Gibbsite(001) 5.2104 5.078 4220 0.939Gibbsite simulations investigated the (001) surface of gibbsite using both six-site and TIP4P/Ice water. Most simulations were carried out at 240 K for bothwater models. Hydroxyl groups were treated in the same way as the free-H simu-lations on kaolinite described in Chapter 5. Hydroxyl hydrogen atoms were ableto pivot around their fixed oxygen atom. We investigated the effect of scaling thelattice parameters of gibbsite to zero the mismatch between gibbsite and kaoliniteor ice. Scalings are summarized in Table 6.2.6.3 ResultsIn all of our simulations we did not observe any ice nucleation from any of thepotassium feldspar surfaces. To understand why no ice nucleated we turn to Figure6.2, where density profiles for all six surfaces we tested are shown for both thesix-site (upper profiles), and TIP4P/Ice (lower profiles), water models at 240 K.In short, these density profiles do not suggest that these surfaces are capable of124Table 6.2: Scalings of the gibbsite (001) surfaces explored. For each scal-ing the lattice parameters are given, along with the match to the latticeconstants of ice Ih.Scaling a×b (A˚) MismatchIha[2:1] Ihc[3:2]Gibbsite(001) 8.684×5.078 3.1% 4.0%Ice-scaled gibbsite(001) 8.96×4.873 0.0% 0.0%Kaolinite-scaled gibbsite(001) 8.9419×5.1535 0.2% 5.7%nucleating ice. In previous studies on silver iodide and kaolinite (Chapters 3 & 5),we observed density profiles in the first layer of water which resembled ice profileseven before ice nucleated. The most common feature, a doublet oxygen peak fromthe formation of a bilayer of ice, with hydrogen peaks between and on either sideof the oxygen doublet, cannot be seen in the profiles shown in Figure 6.2. Thus,we conclude that it is unlikely that these surfaces, at least as they are currentlymodelled, will nucleate ice. A recent combination of DFT and MD simulationsof one to two layers of water on feldspar also suggest that the first layer of wateron feldspar does not resemble ice,[145] although they did find that a second layerof water does become icelike. Our results agree that the first layer of water on thesurface is not icelike, but we do not find icelike ordering in subsequent layers ofwater.The lack of ice nucleation ability of these surfaces conflicts with experimentalevidence suggesting potassium feldspar is an effective IN.[46, 52] We suggest acouple of reasons for this discrepancy:• The ice nucleation ability of feldspar is due to some other surface, or tosome sort of surface defect, or other surface feature that we did not con-12501002003004001.0 1.2 1.4 1.6Z (nm)density (particles/nm3 )(001)01002003004001.1 1.3 1.5 1.7Z (nm)density (particles/nm3 )(001)+O01002003004001.1 1.3 1.5 1.7Z (nm)density (particles/nm3 )(001)+OH02004006001.6 1.8 2.0 2.2Z (nm)density (particles/nm3 )(010)01002003004001.9 2.1 2.3 2.5Z (nm)density (particles/nm3 )(010)+O01002003001.9 2.1 2.3 2.5Z (nm)density (particles/nm3 )(010)+OHFigure 6.2: Density profiles of liquid water near feldspar at 240 K. Wateroxygen and hydrogen profiles are red and black, respectively. Hydro-gen profiles of feldspar hydroxyl groups are green. For each panel,the upper profiles are for six-site water and the lower profiles are forTIP4P/Ice water. Note the density scale for six-site water is shifted 200particles/nm3 for clarity. Above each panel the surface in question isindicated.sider. For instance, it has been suggested that the potassium ions play a rolein ordering water for ice nucleation.[146] In our simulations, potassium ionswere fixed, but giving some freedom to the potassium ions may be neces-sary. Although solvating ions in simulations poses its own problems.[147]• Ice nucleation by feldspar is not caused by the feldspar itself, but by weath-126ering products of feldspar, possibly only present in small amounts on thesurface.The second reason suggested above offers some interesting possibilities. In ex-perimental studies of feldspar, the sample may have been impure.[53, 146, 148]Samples sometimes contain a mixture of several clays and minerals, some beingweathering products of feldspar. Since feldspar cleaves cleanly,[143] it could pro-vide a support to produce smooth surface patches of weathering products that arelarge enough to nucleate ice.Kaolinite is a weathering product of feldspar[58] which we have already stud-ied. Another weathering product of feldspar is gibbsite,[120] which has a surfacethat is structurally very similar to the Al-surface of kaolinite.Our simulations of gibbsite did not produce ice directly. This is shown in Fig-ure 6.3A, where the amount of ice in several simulations is plotted as a function oftime. In an attempt to make the surface more favourable to ice, we considered twoadditional scaled gibbsite surfaces, where the lattice constants have been changedto match either the prism face of ice Ih or kaolinite. Kaolinite was selected be-cause we know that kaolinite nucleates ice, thus we were interested in seeing ifthe gibbsite surface scaled to kaolinite could also nucleate ice.Only the kaolinite-scaled gibbsite simulations nucleated ice under free-H con-ditions (Figure 6.3A). Similar to kaolinite, ice Ih was the dominant form of theice that nucleated and grew. Figure 6.4 shows ice on kaolinite-scaled gibbsite ob-tained with the TIP4P/Ice water model. Ice Ih molecules have their oxygen atomscoloured blue, and the simulation cell is oriented to look down the c-axis of ice127Figure 6.3: Plots of the number of ice molecules for several gibbsite simu-lations. Panel A is for free-H simulations at 240 K with both the six-site (dotted lines), and TIP4P/Ice (solid lines) results included. Redindicates simulations on gibbsite, blue is for ice-scaled gibbsite, andgreen is for kaolinite-scaled gibbsite. Most lines are near the bottom inpanel A because no ice nucleated, making it difficult distinguish them.Panel B is for several ice-prepared simulations. Dotted lines indicatekaolinite-scaled gibbsite and solid lines indicate normal gibbsite. Bluelines are for TIP4P/Ice and red is for six-site water.Ih. Only the TIP4P/Ice model nucleated ice in kaolinite-scaled simulations whenthe water models and hydroxyl groups were coupled to the same temperature. Icenucleation with the six-site model was observed in several simulations where thewater models were regulated to 240 K and the hydroxyl groups regulated to 100K. Reducing the temperature of the hydroxyl groups slows down their dynamics,and apparently helps them to adopt a configuration favourable to ice. This resultfurther shows the role atomistic dynamics plays in simulations of ice nucleation.Kaolinite simulations in Chapter 5 and other simulation work on silver iodide,[83]show that too much movement in a surface can impede ice nucleation on simula-tion timescales.128Figure 6.4: Ice growing on the kaolinite-scaled gibbsite(001) surface withthe TIP4P/Ice water model at 240 K. Oxygen, hydrogen, and alu-minium atoms are red, black, and pink, respectively. Oxygen atomsof ice Ih determined by the CHILL algorithm are blue. The simulationcell is oriented to look down the c-axis of ice Ih.The density profiles in Figure 6.5 provide some explanation for the kaolinite-scaled surface’s ice nucleation ability. Density profiles of ice Ih nucleating onthe prism face of kaolinite (Figures 5.6D and 5.7D) show water oxygen atomsforming a doublet pattern from each bilayer of ice. Each peak in the doublet isabout 1 A˚ apart. The oxygen peak nearest to the surface also shows some splitting.Hydrogen atom density profiles show a large peak in the middle of the oxygendoublet and a small peaks on either side of the doublet. Profiles for each gibbsitesurface investigated (scaled and unscaled), shown in Figure 6.5, are very similar.The most significant difference is in the kaolinite-scaled profiles, where the first129Figure 6.5: Density profiles of liquid water near gibbsite at 240 K. Wateroxygen, water hydrogen, and hydroxyl hydrogen profiles are in red,black, and green respectively. The upper profiles are from six-site wa-ter and the bottom profiles are from TIP4P/Ice water. Note the densityscale for six-site water has been shifted 200 particles/nm3 for clarity.From left to right the plots are for gibbsite, ice-scaled gibbsite, andkaolinite-scaled gibbsite. Zero on the Z-axis coincides with the centreof the surface’s hydroxyl oxygen atoms.oxygen peak at Z = 0.25 nm shows some spitting, and the prominence of thesecond peak of the oxygen doublet at Z ≈ 0.35−0.4 nm. TIP4P/Ice water profilesfor gibbsite and ice-scaled gibbsite (lower profiles in Figure 6.5) are similar to thekaolinite-scaled profiles, but the kaolinite-scaled profiles are more pronounced.From these profiles, it is clear that kaolinite-scaled gibbsite is the most effectivesurface at inducing an icelike configuration in water near its surface.Minor variations in the behaviour of hydroxyl hydrogen atoms may affect icenucleation on simulation time scales. The hydroxyl groups of the first layer ofgibbsite tend to point in one of three directions: toward water and away from thegibbsite slab, horizontally in the plane of oxygen atoms, or away from the wa-ter and into the gibbsite slab. The gibbsite and ice-scaled density profiles shown130in Figure 6.5 demonstrate these three options. In these profiles the hydroxyl hy-drogen density (green curve) has three peaks. From left to right, the first peakrepresents hydroxyl groups pointing into the slab. The middle peak, at zero on thex-axis, represents hydroxyl groups lying horizontally in the plane of the hydroxyloxygen atoms. The third peak shows hydroxyl groups pointing toward water andaway from the slab. Kaolinite-scaled profiles show a much smaller middle peakthan either gibbsite or ice-scaled gibbsite. This suggests that the behaviour of thehydroxyl groups differs for kaolinite-scaled gibbsite compared to other scalingsand this may help ice nucleate on gibbsite. It is remarkable that small differ-ences in scaling could create a significant change in the dynamics of the hydroxylgroups.Figure 6.3B shows the number of ice molecules as a function of time for sev-eral ice-prepared simulations. In these simulations hydroxyl groups are held fixed(including hydroxyl hydrogen atoms) to coordinates obtained after ice had nucle-ated in the kaolinite-scaled simulations. Both gibbsite and kaolinite-scaled gibb-site simulations are shown, and both scalings nucleate ice. Gibbsite ice-preparedsurfaces were created by scaling the coordinates from kaolinite-scaled simula-tions, after ice had nucleated, back to gibbsite’s lattice parameters. This shows forthe rigid surfaces that there is a configuration of gibbsite’s hydroxyl groups ca-pable of nucleating ice. Although not seen in this study, we believe it is possiblefor gibbsite to nucleate ice without resorting to an ice-prepared surface. Observ-ing ice nucleation in direct simulations of gibbsite may be possible, but the timescales required may be quite large and beyond what we are capable of at this time.131Figure 6.6: Density profiles of TIP4P/Ice water (panel A), and six-site water(panel B), near ice-prepared surfaces at 240 K. Upper profiles (note thedensity scale shifted by 200 particles/nm3 for clarity) are for gibbsite,and lower profiles are for kaolinite-scaled gibbsite. Red and blackrepresent oxygen and hydrogen density profiles, respectively. Zero onthe Z-axis coincides with the centre of the surface’s hydroxyl oxygenatoms.Although many IN have been studied experimentally, we are unaware of any ex-perimental studies of ice nucleation by gibbsite to compare with our results.An explanation of the ice-prepared surface’s ice nucleation ability can be seenin the density profiles in Figure 6.6 of water near ice-prepared surfaces. Profilesfor both water models near gibbsite and kaolinite-scalings are shown. Densityprofiles for TIP4P/Ice (panel A) show the oxygen doublets forming for both scal-ings and clear splitting in the first peak, demonstrating that the surface is effectiveat inducing an icelike configuration on water near the surface. Six-site water inpanel B does not show as clear a splitting in the first peak as in panel A. However,the surfaces in panel B are capable of nucleating ice, so the surface induced orderis sufficient to nucleate ice.1326.4 ConclusionSimulations were conducted on feldspar to determine its mechanism of ice nucle-ation. The feldspar surfaces tested did not nucleate ice and the water density pro-files suggest that these surfaces are unlikely to nucleate ice, at least as we modelthem. We hypothesized that the ice nucleation ability of feldspar may be due, inpart, to weathering products on its surface, such as kaolinite and/or gibbsite. The(001) surface of gibbsite is very similar to the (001) surface of kaolinite which wehave shown to nucleate ice. We were unable to obtain ice nucleation in direct sim-ulations of gibbsite, however, by scaling gibbsite to match the lattice parametersof kaolinite we observed ice nucleation in several simulations. Ice nucleating onthe kaolinite-scaled gibbsite ordered the surface hydroxyl groups into a configu-ration compatible with ice. Simulations with this ice-prepared surface held fixednucleated ice with both kaolinite-scaled gibbsite and with normal gibbsite latticeparameters. These ice-prepared simulations show that there exists a configurationof gibbsite which can nucleate ice.133Chapter 7Attempted Ice Nucleation by aProtein“What is the secret of life?” I asked.“I forget,” said Sandra.“Protein,” the bartender declared. “They found out something aboutprotein.”“Yeah,” said Sandra, “thats it.”— Kurt Vonnegut, Cat’s Cradle7.1 SummaryBiological aerosols are thought to represent a significant source of IN in the atmo-sphere. This project’s goal was to simulate ice nucleation by a protein. Ultimatelythe goal of this project was not achieved and possible reasons why are discussed.This chapter is included in this thesis to document what we have done, and thethinking behind it, as a resource to anyone who may attempt this project in thefuture.1347.2 StrategyOne of the most experimentally studied ice nucleation proteins (INP) comes fromthe bacteria Pseudomonas syringae. The ice nucleation efficiency of this bacteriarivals that of silver iodide because it can nucleate ice at temperatures as warm as -3◦C.[59] Isolation of the gene responsible for the INP revealed that most of the se-quence of amino acids repeats every sixteen residues.[60] The weblogo plot[149]in Figure 7.1 shows the probabilities of amino acids for residues 208-959 out of1200 total residues in the INP. A weblogo plot displays repetition in a sequenceof amino acids. Here, the probabilities of an amino acid appearing every sixteenresidues is shown. Figure 7.1 shows that the amino acids at positions 1-6, 8, 9,13, and 15 of every sixteen residue unit of the P. syringae INP are very consistent,with lesser amounts of consistency at other positions. This consistency shows therepetitive nature of this portion of the protein. This repeating section may takethe shape of a β -helix,[150] which is a helix forming back-to-back β -sheets. Theamino acid sequence at the ends of the protein deviate from this repeating patten.Clusters of the P. syringae INP are responsible for ice nucleation, as many as50-60 proteins are required to nucleate ice at -3◦C, and only a single protein is re-quired at -13◦C.[61] A recent combination of simulation and experimental workhas shown that the P. syringae INP can order water on its surface, and that theprotein may have some heat dissipation ability which helps ice nucleation.[151]Simulations of the INP from Pseudomonas borialis, which has a very similar se-quence to the P. syringae INP, show that the P. borialis INP can dimerize, creating135positonWebLogo 3.40.00.51.0probabilitySTATGYG5STASGQTSA10SERQGDVQHFGAYSETQHSKRNAEGDTSWVMFIDATNS15ILVITFigure 7.1: Weblogo[149] plot of the probabilities of residues in the P.syringae INP for residues 208-959. Position is the residue lo-cation in a repeated sixteen residue segment. Letters representamino acids, and the size of the letter shows the probability atwhich that amino acid appears at that location in the repeated seg-ment. Blue indicates hydrophillic residues, green is neutral, andblack is hydrophobic. Residues appearing here are: A=Alanine,D=Aspartic acid, E=Glutamic acid, F=Phenylalanine, G=Glycine,H=Histidine, I=Isoleucine, K=Lysine, L=Leucine, M=Methionine,N=Asparagine, Q=Glutamine, R=Arginine, S=Serine, T=Threonine,V=Valine, W=Tryptophan, Y=Tyrosine.a larger flat surface from the combined β -sheets.[62] Water on the surface of theβ -sheet localizes into a pattern which matches ice.[62]Because there is no crystal structure of the P. syringae INP, our strategy wasto find a protein with a β -helix structure, ideally similar in size to the P. syringaeINP and with a repetitive sequence compatible with ice, and mutate the proteinto create simulations with several proteins aligned to form a raft. These simula-tions were intended to mimic the slab simulations in previous chapters which weresuccessful in nucleating ice.1367.3 Methods and ResultsA search of the RCSB Protein Data Bank[116] revealed several candidate proteins,and we selected the 3ULT[121] protein, which is an anti-freeze protein (AFP)from the rye grass Lolium perenne. Many plants and animals in cold climates havedeveloped AFPs.[152] AFPs often bear similarities to INPs because AFPs tend toconsist of repeating helices, although usually shorter than INPs. Experimentalstudies have shown that the L. perenne AFP binds well to the basal and prismfaces of Ih.[121]The original 3ULT protein is a helix (Figure 7.2A), containing a fourteenresidue repeating pattern, forming two β -sheets. One sheet has three loops con-taining an extra amino acid creating an outward protrusion. Experiments haveshown this protruding side has inferior ice binding abilities compared to the othersmooth side.[153] Our mutations of this protein removed the extra amino acids,and some residues were changed to make each β -sheet more uniform by havingeach β -sheet contain two hydrophilic residues whose side chains pointed awayfrom the protein. Several of the residues changed were hydrophobic, and theiralteration resulted in a more uniform hydrophillic surface. This regular threonine-valine-serine sequence is highlighted in Figure 7.2B. The side chains of threonineand serine extend away from the protein and the valine side chain resides insidethe helix. Two of these mutated proteins were then connected to make the finalprotein longer.In attempting ice nucleation, we used a configuration of three mutated proteins137­PNTISGSNNTVRSGSKNVLAGNDNTVISGDNNSVSGSNNTVVSGNDNTVTGSNHVV­SGTNHIVTDNNNNV­SGNDNNVSGSFHTV­SGGHNTVSGSNNTV­SGSNHVVSGSNKVVTD­­­­­­­SNNTVSGSKNTVSGNDNTVSGDNNTVSGSNNTVSGNDNTVSGSNHTVSGTNHTVSDNNNTVSGNDNTVSGSFHTVSGGHNTVSGSNNTVSGSKNTVSGNDNTVSGDNNTVSGSNNTVSGNDNTVSGSNHTVSGTNHTVSDNNNTVSGNDNTVSGSFHTVSGGHNTVSGA                        BFigure 7.2: Sequences and images of the proteins used in this study. PanelA displays the sequence (top panel) of the original AFP protein fromL. perenne, and an image (bottom panel) of the protein. Residues inthe sequence shown in red were removed during mutation, and theresidues in blue were changed. Panel B displays the sequence andimage of the mutated protein used in the simulations. Residues shownin blue were changed from the original AFP. In the image in Panel Bthreonine residues are red, valine residues are blue, and serine residuesare orange.138such that they produced a repeating surface in the X/Y plane of the simulation cell.Infinite slab simulations were successful at nucleating ice in preceding chapters,so we reasoned a similar geometry could be successful with proteins. Figure 7.3shows this configuration, which we refer to as a raft. The proteins are lined upparallel to each other, and will connect seamlessly to the periodic images in theX and Y directions. Test simulations we conducted showed that the 3ULT proteinwould form dimers when aligned parallel to each other, so we feel that this sortof configuration is reasonable. Alternative alignments, such as anti-parallel, didnot remain together or create extended flat surfaces like the parallel arrangement.A simulation run at 1 bar pressure and 240 K for six-site water, and at 230 K forTIP4P/Ice water, did not show any sign of ice nucleation for either water model.The amber03 force field[124] was used to parametrize the protein.Experience with kaolinite, Chapter 5, suggested that a surface may be ableto adopt certain configurations which allow it to nucleate ice quite readily. Weattempted to induce such an icelike configuration on our protein by simulatingthe protein in contact with ice generated by an AgI disk from Chapter 4 withsix-site water. Disk E, was used for this because it nucleated ice at the warmesttemperatures. Only one layer of the disk was used and the charges on AgI were setto zero so there would be no electric field created by the disk. Two orientations ofthe protein were considered, one where the protein was oriented perpendicular tothe surface of the disk, and another where the protein ran approximately parallelto the disk (Figure 7.4). Each simulations was run under NPT conditions at 1 barpressure, and at a temperature of 260 K.139ABZ    XY    XFigure 7.3: Two views of the raft of proteins. Panel A shows a top viewand panel B shows an end-on view. The proteins are situated in thesimulation cell to produce an infinitely periodic surface of proteins.Unfortunately, the ice nucleated on the surface of the AgI disk did not com-pletely reach the protein. A thin layer of liquid water continued to surround theprotein, thus the ice did not induce an icelike configuration on the protein. A sim-ulation using these “ice prepared” proteins in the raft configuration was createdfor each water model and run at 220 K. The protein was held fixed in these sim-ulations so as not to destroy the configuration. The simulations did not yield anyice.A closer look at the behaviour of water around the proteins during the AgI disksimulations shows how the protein might nucleate ice. In Figure 7.5 the regions140A BFigure 7.4: Two orientations of the protein used in an attempt to “force”the protein into an icelike configuration. The protein (teal) is orientedeither approximately parallel to the AgI disk (green and silver) as inpanel A, or perpendicular to the disk as in panel B. Water moleculesare not shown.occupied by water oxygen atoms around the protein are shown in pink. Threonine,valine, and serine residues are also shown. In both cases, most clearly observedin panel B, the threonine and serine residues produce a region without water. Oneither side of the threonine residues are channels with water extending towardthe protein. Creation of these two channels of water likely restricts some of themovement of the water, which may help water molecules to settle into an icelikeconfiguration.141A                      BFigure 7.5: Horizontal (Panel A) and vertical (Panel B) “ice-prepared” pro-teins. Colours of the proteins are the same as in Figure 7.2, and thestructure of threonine, valine, and serine residues is displayed. Pinkshapes show the regions occupied by water oxygen atoms.7.4 ConclusionsIn the end, this project failed in its goal to obtain a model protein-like surface thatfunctions as an IN. One can think of several possible reasons for this. For exam-ple, the secondary structure of the protein may not resemble the P. syringae INP,an ineffective sequence was used, or the simulation geometry did not accuratelyrepresent what happens with real proteins. It would be useful for more experi-mental evidence to be collected before this project is attempted again. It would beparticularly useful to know the crystal structure of the P. syringae INP.142Chapter 8Conclusion“If I were a younger man, I would write a history of human stupidity;and I would climb to the top of Mount McCabe and lie down on myback with my history for a pillow; and I would take from the groundsome of the blue-white poison that makes statues of men; and I wouldmake a statue of myself, lying on my back, grinning horribly, andthumbing my nose at You Know Who.”— Kurt Vonnegut, Cat’s Cradle8.1 Summary of ThesisWater droplets do not freeze homogeneously above -38◦C without the aid of anIN.[1] Many different IN have been discovered, each with different temperatureranges over which they are active.[1] What is lacking is an understanding of themicroscopic mechanisms by which IN nucleate ice. This thesis attempts to dis-cover the microscopic mechanisms of ice nucleation by several IN using molecu-lar dynamics simulations.The first IN investigated was silver iodide. In Chapter 4, we studied seven143surfaces from two phases of silver iodide. From β -AgI we tested the silver andiodide exposed (0001) surface and the Ag+I exposed (101¯0) surface. From γ-AgIwe tested the silver and iodide exposed (001) and (111) surfaces. Of the surfaceswe looked at, we found that only the silver exposed surfaces nucleated ice. Uponinspection of these surfaces via images and density profiles, we found silver ex-posed surfaces organized water into an icelike configuration, even at temperatureswarmer than where ice nucleation occurred. Favourable configurations includedchair conformed rings and face-centred cubic arrangements with water moleculesoriented in a similar fashion to bulk ice. Other surfaces investigated did not or-ganize water into an icelike configuration, and did not nucleate ice. Non-icelikeconfigurations included planar rings which are not found in bulk ice.Chapter 4 explored finite silver iodide disks and plates. Disks and plates couldbe representative of small patches of a smooth surface, which may exist on largersilver iodide particles. Five disks and plates with radii between 1.15 to 2.99 nmwere considered. Larger disks were found to nucleate ice at warmer temperaturesthan smaller plates, approaching asymptotically the temperature of nucleation ob-served on infinite silver iodide surfaces. This trend of larger disks nucleating ice atwarmer temperatures is consistent with previous simulations and theoretical pre-dictions. Indeed, there is a very good match with the GT equation. The geometryof the disks and plates did not affect the results. Ice clusters on the disks oftengrew and shrank a few times before reaching some critical size followed by rapidgrowth of bulk ice. Estimates of the size of the critical clusters in our simula-tions is consistent with trends from previous theoretical and experimental studies.144The effect of the disks is to provide an initial cluster of ice which has difficultlymelting, and thus the critical cluster size is easier to reach.Kaolinite is a common clay which has long been known as an effective IN.[1]Simulations described in Chapter 5 focused on both sides of the (001) surface. TheAl-surface is comprised of hydroxyl groups. We simulated the hydroxyl groupswith varying amounts of freedom of motion which we termed rigid-kaolinite,free-H, and free-OH. Rigid-kaolinite held kaolinite’s hydroxyl groups fixed totheir crystallographic coordinates, free-H allowed the hydroxyl hydrogen atomsto move around their oxygen atoms, and free-OH also allowed the hydroxyl oxy-gen atoms to vibrate around their crystallographic coordinates. Ice only nucleatedin free-H simulations, and required 100-300 ns before nucleating, and only ice Ihformed. Some freedom of movement in the kaolinite surface was found necessaryfor kaolinite to adapt into a configuration capable of nucleating ice. Too muchfreedom (i.e. movement) prevented the surface from reaching an icelike config-uration on our simulation time scales. Once a surface nucleated ice, the result-ing configuration, which had an identifiable pattern, of the hydroxyl groups wascompatible with ice and proved to be an efficient surface for nucleating ice. Con-figurations of hydroxyl groups taken before ice had formed, were not observed tonucleate ice.The Si-surface of kaolinite was also found to nucleate ice; however, this sur-face did so by organizing the first layer of water into a configuration which doesnot resemble any phase of ice we are aware of. The first layer of ice contained arepeating pattern of five, six, and eight member rings which grew into a regular145mixture of Ih and Ic. This result suggests that water is capable of forming novelstructures on a surface, which can facilitate ice nucleation.Feldspar has recently been discovered as a very effective IN.[46, 52] We stud-ied (001) and (010) surfaces of potassium feldspar in Chapter 6. Three variationsof each surface were tested: bare, containing only silicon, aluminium, potassium,and bridging oxygen atoms; +O, which had additional oxygen atoms; and +OH,which had additional hydroxyl groups. None of these surfaces nucleated ice, anddensity profiles from these simulations suggest that these surfaces will not be ableto nucleate ice. We then turned to gibbsite, a weathering product of feldspar.[54–57] Gibbsite shares many surface features with kaolinite. No ice nucleation wasobserved on gibbsite, however we were able to nucleate ice by re-scaling gibb-site’s lattice parameters to match kaolinite. Scaled simulations produced surfaceswhich were complementary to ice, which we then scaled back to gibbsite and usedto nucleate ice. This results shows there is a configuration of the gibbsite surfacewhich can nucleate ice on simulation timescales.Finally in Chapter 7 we attempted to nucleate ice using a protein. No crys-tal structures of INPs exist so we used an anti-freeze protein structure in oursimulations.[121] The protein was mutated and arranged to produce an infiniteraft of proteins. These simulations did not nucleate ice, even after attempting toinduce an icelike configuration on the protein.The work in this thesis sheds light on possible mechanisms of ice nucleationfor silver iodide, kaolinite, and gibbsite. Ice nucleation mechanisms on feldsparand proteins are still elusive as are the mechanisms by which other IN operate.1468.2 Future WorkThere is a lack of experimental evidence on what the surfaces of IN in the atmo-sphere or laboratory look like at the atomistic level. Therefore it would be benefi-cial for more experimental data to be available on the surfaces present for variousIN. Simulations could then be configured to be consistent with this information.Surface structure information from electron microscopy, atomic force microscopy,scanning tunnelling microscopy, and x-ray techniques would be valuable.The majority of the studies in this thesis examined infinite flat surfaces. On theother hand, aerosol particles might not have very large flat patches on them. Verylarge patches are not necessary for ice nucleation, but we do not know much aboutthe sizes or types of surfaces which do exist. Future work should examine surfacedefects and other structures which can occur on an IN’s surface. Particularly forprotein simulations, the availability of the crystal structure of one or more INPswould be useful.147Bibliography[1] Pruppacher, H. R.; Klett, J. D. Microphysics of Clouds and Precipitation.Springer, New York, 2010. → pages 2, 5, 6, 7, 8, 9, 44, 106, 143, 145[2] Sanz, E.; Vega, C.; Espinosa, J. R.; Caballero-Bernal, R.; Abascal, J. L. F.;Valeriani, C. Homogeneous Ice Nucleation at Moderate Supercoolingfrom Molecular Simulation. J. Am. Chem. Soc. 2013. 135, 15008–15017.→ pages 2, 12, 74, 85, 86, 87[3] Seinfeld, J. H.; Pandis, S. N. Atmospheric Chemistry and Physics: FromAir Pollution to Climate Change, 2nd ed. Wiley-Interscience, Hoboken,N.J., 2006. → pages 2, 5, 6, 7, 68[4] Hegg, D. A.; Baker, M. B. Nucleation in the Atmosphere. Rep. Prog.Phys. 2012. 72, 056801. → pages 2, 4[5] Espinosa, J. R.; Vega, C.; Valeriani, C.; Sanz, E. Seeding Approach toCrystal Nucleation. J. Chem. Phys. 2016. 144, 034501. → pages 2[6] Molinero, V.; Moore, E. B. Water Modeled As an Intermediate ElementBetween Carbon and Silicon. J. Phys. Chem. B 2009. 113, 4008–4016. →pages 3, 9, 12, 13, 37, 74, 85[7] Ickes, L.; Welti, A.; Hoose, C.; Lohmann, U. Classical Nucleation Theoryof Homogeneous Freezing of Water: Thermodynamic and KineticParameters. Phys. Chem. Chem. Phys. 2015. 17, 5514–5537. → pages 4[8] Fletcher, N. H. Size Effect in Heterogeneous Nucleation. J. Chem. Phys.1958. 29, 572–576. → pages 5[9] Fletcher, N. H. Nucleation by Crystalline Particles. J. Chem. Phys. 1963.38, 237–240. → pages 5148[10] Fletcher, N. H. On Ice-Crystal Production by Aerosol Particles. J.Meteorol. 1959. 16, 173–180. → pages 5[11] Stocker, T. F.; Qin, D.; Plattner, G.-K.; Tignor, M. M. B.; Allen, S. K.;Boschung, J.; Nauels, A.; Xia, Y.; Bex, V.; Midgley, P. M. IPCC, 2013:Climate Change 2013: The Physical Science Basis. Contribution ofWorking Group I to the Fifth Assessment Report of the IntergovernmentalPanel on Climate Change. Cambridge University Press, Cambridge, UKand New York, NY, 2013. → pages 6, 10[12] Lau, K. M.; Wu, H. T. Warm Rain Processes Over Tropical Oceans andClimate Implications. Geophys Res Lett 2003. 30, 2290–2295. → pages 6[13] Murray, B. J.; O’Sullivan, D.; Atkinson, J. D.; Webb, M. E. Ice Nucleationby Particles Immersed in Supercooled Cloud Droplets. Chem Soc Rev2012. 41, 6519–6554. → pages 6, 10[14] Ervens, B.; Feingold, G. On the Representation of Immersion andCondensation Freezing in Cloud Models Using Different NucleationSchemes. Atmos Chem Phys 2012. 12, 5807–5826. → pages 7[15] Turnbull, D.; Vonnegut, B. Nucleation Catalysis. Ind. Eng. Chem. 1958.44, 1292–1298. → pages 8[16] Pinti, V.; Marcolli, C.; Zobrist, B.; Hoyle, C. R.; Peter, T. Ice NucleationEfficiency of Clay Minerals in the Immersion Mode. Atmos. Chem. Phys.2012. 12, 5859–5878. → pages 8, 10[17] Bartels-Rausch, T.; Bergeron, V.; Cartwright, J. H. E.; Escribano, R.;Finney, J. L.; Grothe, H.; Gutie´rrez, P. J.; Haapala, J.; Kuhs, W. F.;Pettersson, J. B. C.; Price, S. D.; Sainz-Dı´az, C. I.; Stokes, D. J.;Strazzula, G.; Thomson, E. S.; Trinks, H.; Uras-Aytemiz, N. IceStructures, Patterns, and Processes: A View Across the Icefields. Rev.Mod. Phys. 2012. 84, 885–944. → pages 8[18] Malenkov, G. Liquid Water and Ices: Understanding the Structure andPhysical Properties. J. Phys: Condens. Matter. 2009. 21, 283101–1–35.→ pages 8, 9, 107149[19] Russo, J.; Romano, F.; Tanaka, H. New Metastable Form of Ice and itsRole in the Homogeneous Crystallization of Water. Nat. Matter. 2014. 13,733–739. → pages 8[20] Mayer, E.; Hallbrucker, A. Cubic Ice from Liquid Water. Nature 1987.325, 601–602. → pages 8[21] Murray, B. J.; Knopf, D. A.; Bertram, A. K. The Formation of Cubic IceUnder Conditions Relevant to Earths Atmosphere. Nature 2005. 434,202–205. → pages 8[22] Carignano, M. A. Formation of Stacking Faults During Ice Growth onHexagonal and Cubic Substrates. J. Phys. Chem. C 2007. 111, 501–504.→ pages 8, 9[23] Johari, G. P. On the Coexistence of Cubic and Hexagonal Ice Between160 and 240 K. Philos. Mag. 1998. 78, 375–383. → pages 9, 45, 118[24] Malkin, T. L.; Murray, B. J.; Brukhno, A. V.; Anwar, J.; Salzmann, C. G.Structure of Ice Crystallized from Supercooled Water. P. Natl. Acad. Sci.USA 2012. 109, 1041–1045. → pages 9, 45, 116[25] Kuhs, W. F.; Sippel, C.; Falenty, A.; Hansen, T. C. Extent and Relevanceof Stacking Disorder in “Ice Ic”. P. Natl. Acad. Sci. USA 2012. 109,21259–21264. → pages 9, 45[26] Malkin, T. L.; Murray, B. J.; Salzmann, C. G.; ; Molinero, V.; Pickering,S. J.; Whale, T. F. Stacking Disorder in Ice I. Phys. Chem. Chem. Phys.2015. 17, 60–76. → pages 9, 72, 116[27] Bruninjes, R. T. A Review of Cloud Seeding Experiments to EnhancePrecipitation and Some New Prospects. B Am Meteorol Soc 1999. 80. →pages 9[28] Breed, D.; Rasmussen, R.; Weeks, C.; Boe, B.; Deshler, T. EvaluatingWinter Orogrophic Cloud Seeding: Design of the Wyoming WeatherModification Pilot Project (WWMPP). J. Appl. Meteorol. Clim. 2014. 53,282–299. → pages 9[29] Vonnegut, B. The Nucleation of Ice Formation by Silver Iodide. J ApplPhys 1947. 18, 593–595. → pages 9, 44150[30] Zettlemoyer, A. C.; Tcheurekdjian, N.; Chessick, J. J. Surface Propertiesof Silver Iodide. Nature 1961. 192, 653. → pages 9[31] Finnegan, W. G.; Chai, S. K. A New Hypothesis for the Mechanism of IceNucleation on Wetted AgI and AgI-AgCl Particulate Aerosols. J AtmosSci 2003. 60, 1723–1731. → pages 9[32] Burley, G. Polymorphism of Silver Iodide. Am. Mineral. 1963. 48,1266–1276. → pages 9[33] Burley, G. Structure of Hexagonal Silver Iodide. J Chem Phys 1963. 38,2807–2812. → pages 9, 30, 64[34] Fukata, N.; Paik, Y. Water Adsorption and Ice Nucleation on Silver IodideSurfaces. J. App. Phys. 1973. 44, 1092–1100. → pages 10, 48[35] Hale, B. N.; Kiefer, J. Studies of H2O on β -AgI Surfaces: An EffectivePair Potential Model. J Chem Phys 1980. 73, 923–933. → pages 30, 172,173[36] Hale, B. N.; Kiefer, J.; Terrazas, S.; Ward, R. C. Theoretical Studies ofWater Adsorbed on Silver Iodide. J. Phys. Chem. 1980. 84, 1473–1479.→ pages[37] Shevkunov, S. V. Adsorption of Water Vapor on the AgI Surface: AComputer Experiment. Russ. J. Gen. Chem. 2005. 75, 1632–1642. →pages[38] Shevkunov, S. V. Layer-by-Layer Adsorption of Water Molecules on theSurface of a Silver Iodide Crystal. Russ. J. Phys. Ch. 2006. 80, 769–775.→ pages[39] Shevkunov, S. V. Clusterization of Water Molecules on Crystalline β -AgISurface. Computer Experiment. Colloid J+ 2006. 68, 632–643. → pages[40] Taylor, J. H.; Hale, B. N. Monte Carlo Simulations of Water-Ice Layers ona Model Silver Iodide Substrate: A Comparison with Bulk Ice Systems.Phys. Rev. B. 1993. 47, 9732–9741. → pages 51[41] Volkov, V. B.; Zhogolev, D. A.; Bakhanova, R. A.; Kiselev, V. J. AQuantum-Chemical Study of the Interaction of the Water Molecule with151Silver Iodide Surface Clusters. Chem. Phys. Lett. 1976. 43, 45–48. →pages[42] Ward, R. C.; Holdman, J. M.; Hale, B. N. Monte Carlo Studies of WaterMonolayer Clusters on Substrates: Hexagonal AgI. J. Chem. Phys. 1982.77, 3198–3202. → pages[43] Ward, R. C.; Hale, B. N.; Terrazas, S. A Study of the Critical Cluster Sizefor Water Monolayer Clusters on a Model AgI Basal Substrate. J. Chem.Phys. 1983. 78, 420–423. → pages 10, 48[44] Hoose, C.; Mo¨hler, O. Heterogeneous Ice Nucleation on AtmosphericAerosols: A Review of Results from Laboratory Experiments. Atmos.Chem. Phys 2012. 12, 9817–9854. → pages 10[45] Bartels-Rausch, T. Ten Things We Need to Know About Ice and Snow.Nature 2013. 494, 27–29. → pages 10[46] Atkinson, J. D.; Murray, B. J.; Woodhouse, M. T.; Whale, T. F.; Baustian,K. J.; Carslaw, K. S.; Dobbie, S.; O’Sullivan, D.; Malkin, T. L. TheImportance of Feldspar for Ice Nucleation by Mineral Dust inMixed-Phase Clouds. Nature 2013. 498, 355–358. → pages 10, 121, 125,146[47] Zimmermann, F.; Ebert, M.; Worringen, A.; Schu¨tz, L.; Weinbruch, S.Environmental Scanning Electron Microscopy (ESEM) as a newTechnique to Determine the Ice Nucleation Capability of IdividualAtmospheric Aerosol Particles. Atmos. Environ. 2007. 41, 8219–8227. →pages 10[48] Zimmermann, F.; Weinbruch, S.; Schu¨tz, L.; Hofmann, H.; Ebert, M.;Kandler, K.; Worringen, A. Ice Nucleation Properties of the MostAbundant Mineral Dust Phases. J. Geophys. Res 2008. 113,D23204–1–11. → pages 10[49] Anthony, J. W.; Bideaux, R. A.; Bladh, K. W.; Nichols, M. C. Handbookof Mineralogy. Mineral Data Publishing, Chantilly, VA, 2001. → pages10, 122152[50] S˘olc, R.; Gerzabek, M. H.; Lischka, H.; Tunega, D. Wettability ofKaolinite (001) Surfaces - Molecular Dynamic Study. Geoderma 2011.169, 47–54. → pages 10, 107, 113[51] Tunega, D.; Gerzabek, M. H.; Lischka, H. Ab Initio Molecular DynamicsStudy of a Molecular Water Layer on Octahedral and TetrahedralKaolinite Surfaces. J. Phys. Chem. B 2004. 108, 5930–5936. → pages 10,14, 105[52] Yakobi-Hancock, J. D.; Ladino, L. A.; Abbatt, J. P. D. Feldspar Mineralsas Efficient Deposition Ice Nuclei. Atmos. Chem. Phys. 2013. 13,11175–11185. → pages 10, 121, 125, 146[53] Augustin-Bauditz, S.; Wex, H.; Kanter, S.; Ebert, M.; Niedermeier, D.;Stolz, F.; Prager, A.; Stratmann, F. The Immersion Mode Ice NucleationBehavior of Mineral Dusts: A Comparison of Different Pure and SurfaceModified Dusts. Geophys. Res. Lett. 2014. 41, 7375–7382. → pages 11,127[54] Gardner, L. R. A Chemical Model of the Origin of Gibbsite fromKaolinite. Am. Mineral. 1970. 55, 1380–1389. → pages 11, 146[55] Gardner, L. R. Conditions for Direct Formation of Gibbsite fromK-Feldspar–Further Discussion. Am. Mineral. 1972. 294–300. → pages[56] Lodding, W. Conditions for Direct Formation of Gibbsite fromK-Feldspar–Discussion. Am. Mineral. 1972. 57, 292–294. → pages[57] Kawano, M.; Tomita, K. Amorphous Aluminum Hydroxide Formed at theEarliest Weathering Stages of K-Feldspar. Clay. Clay. Miner. 1996. 44,672–676. → pages 146[58] Wilson, M. J. Weathering of the Primary Rock-Forming Minerals:Processes, Products and Rates. Clay Miner. 2004. 39, 233–266. → pages11, 127[59] Maki, L. R.; Galyan, E. L.; Chang-Chien, M.; Caldwell, D. R. IceNucleation Induced by Pseudomonas Syringae. Appl. Microbiol. 1974.28, 456–459. → pages 11, 135153[60] Green, R. L.; Warren, G. J. Physical and Functional Repitition in aBacterial Ice Nucleation Gene. Nature 1985. 317, 645–648. → pages 11,135[61] Govindarajan, A. G.; Lindow, S. E. Size of Bacterial Ice-Nucleation SitesMeasured in situ by Radiation Inactivation Analysis. Proc. Natl. Acad.Sci. 1988. 85, 1334–1338. → pages 11, 135[62] Garnham, C. P.; Campbell, R. L.; Walker, V. K.; Davies, P. L. NovelDimeric β -Helical Model of an Ice Nucleation Protein with BridgedActive Sites. BMC Struct. Biol. 2011. 11, 36. → pages 11, 136[63] Sosso, G. C.; Chen, J.; Cox, S. J.; Fitzner, M.; Pedevilla, P.; Zen, A.;Michaelides, A. Crystal Nucleation in Liquids: Open Questions andFuture Challenges in Molecular Dynamics Simulations. Chem. Rev. 2016.116, 7078–7116. → pages 12[64] Espinosa, J. R.; Sanz, E.; Valeriani, C.; Vega, C. Homogeneous IceNucleation Evaluated for Several Water Models. J. Chem. Phys. 2014.141, 18C529–1–14. → pages 12, 74, 85, 86[65] Pereyra, R. G.; Carignano, M. A. Ice Nanocolumns: A MolecularDynamics Study. J. Phys. Chem. C 2009. 113, 12699–12705. → pages 74[66] Pereyra, R. G.; Szleifer, I.; Carignano, M. A. Temperature Dependence ofIce Critical Nucleus Size. J. Chem. Phys. 2011. 135, 034508–1–5. →pages 12, 74, 76[67] Matsumoto, M.; Saito, S.; Ohmine, I. Molecular Dynamics Simulation ofthe Ice Nucleation and Growth Process Leading to Water Freezing. Nature2002. 416, 409–413. → pages 12[68] Moore, E. B.; Molinero, V. Structural Transformation in SupercooledWater Controls the Crystallization Rate of Ice. Nature 2011. 479,506–509. → pages 12[69] Stillinger, F. H.; Weber, T. A. Computer Simulation of Local Order inCondensed Phases of Silicon. Phys. Rev. B 1985. 31, 5262–5271. →pages 12154[70] Haji-Akbari, A.; Debenedetti, P. G. Direct Calculation of IceHomogeneous Nucleation Rate for a Molecular Model of Water. Proc.Natl. Acad. Sci. 2015. 112, 10582–10588. → pages 12[71] Brukhno, A. V.; Anwar, J.; Davidchack, R.; Handel, R. Challenges inMolecular Simulation of Homogeneous Ice Nucleation. J. Phys.:Condens. Matter 2008. 20, 494243. → pages 12[72] Abascal, J. L. F.; Vega, C. A General Purpose Model for the CondensedPhases of Water: TIP4P/2005. J. Chem. Phys. 2005. 123, 234505–1–12.→ pages 12, 29, 30, 85[73] Svishchev, I. M.; Kusalik, P. G. Crystallization of Liquid Water in aMolecular Dynamics Simulation. Phys. Rev. Lett. 1994. 73, 975–978. →pages 13, 33[74] Yan, J. Y.; Patey, G. N. Heterogeneous Ice Nucleation Induced by ElectricFields. J. Phys. Chem. Lett. 2011. 2, 2555–2559. → pages 13, 37[75] Yan, J. Y.; Patey, G. N. Molecular Dynamics Simulations of IceNucleation by Electric Fields. J. Phys. Chem. A 2012. 116, 7057–7064.→ pages 41, 42, 67[76] Yan, J. Y.; Patey, G. N. Ice Nucleation by Electric Surface Fields ofVarying Range and Geometry. J. Chem. Phys. 2013. 139, 144501. →pages 13[77] Yan, J. Y.; Overduin, S. D.; Patey, G. N. Understanding Electrofreezing inWater Simulations. J. Chem. Phys. 2014. 141, 074501–1–9. → pages 13,85[78] Braslavsky, I.; Lipson, S. G. Electrofreezing Effect and Nucleation of IceCrystals in Free Growth Experiments. Appl. Phys. Lett. 1998. 72,264–266. → pages 13[79] Stan, C. A.; Tang, S. K. Y.; Bishop, K. J. M.; Whitesides, G. M.Externally Applied Electric Fields Up to 1.6 A˚10 5 V/m Do Not Affect theHomogeneous Nucleation of Ice in Supercooled Water. J. Phys. Chem. B2011. 115, 1089–1097. → pages 13155[80] Aragones, J. L.; MacDowell, L. G.; Siepmann, J. I.; Vega, C. PhaseDiagram of Water Under an Applied Electric Field. Phys. Rev. Lett. 2011.107, 155702. → pages 13[81] Lupi, L.; Hudait, A.; Molinero, V. Heterogeneous Nucleation of Ice onCarbon Surfaces. J. Am. Chem. Soc. 2014. 136, 3156–3164. → pages 13,74[82] Lupi, L.; Molinero, V. Does Hydrophilicity of Carbon Particles ImproveTheir Ice Nucleation Ability? J. Phys. Chem. A 2014. 118, 7330–7337.→ pages 13[83] Fraux, G.; Doye, J. P. K. Note: Heterogeneous Ice Nucleation onSilver-Iodide-Like Surfaces. J. Chem. Phys. 2014. 141, 216101–1–2. →pages 13, 62, 70, 94, 128[84] Cox, S. J.; Kathmann, S. M.; Slater, B.; Michaelides, A. MolecularSimulations of Heterogeneous Ice Nucleation. I. Controlling IceNucleation Through Surface Hydrophilicity. J. Chem. Phys. 2015. 142,184704–1–5. → pages 13[85] Cox, S. J.; Kathmann, S. M.; Slater, B.; Michaelides, A. MolecularSimulations of Heterogeneous Ice Nucleation. II. Peeling Back the Layers.J. Chem. Phys. 2015. 142, 184705–1–8. → pages 13[86] Cox, S. J.; Raza, Z.; Kathmann, S. M.; Slater, B.; Michaelides, A. TheMicroscopic Features of Heterogenous Ice Nucleation May Affect theMacroscopic Morphology of Atmospheric Ice Crystals. Faraday. Discuss.2013. 167, 389–403. → pages 14, 98, 107[87] Croteau, T.; Bertram, A. K.; Patey, G. N. Adsorption and Structure ofWater on Kaolinite Surfaces: Possible Insight into Ice Nucleation fromGrand Canonical Monte Carlo Calculations. J Phys Chem A Lett 2008.112, 10708–10712. → pages 14, 98, 105[88] Croteau, T.; Patey, A. K. B. G. N. Simulation of Water Adsorption onKaolinite Under Atmospheric Conditions. J Phys Chem 2009. 113,7826–7833. → pages 14, 33, 67, 98156[89] Croteau, T.; Bertram, A. K.; Patey, G. N. Water Adsorption on KaoliniteSurfaces Containing Trenches. J. Phys. Chem. A 2010. 114, 2171–2178.→ pages 14[90] Croteau, T.; Bertram, A. K.; Patey, G. N. Observations of High-DensityFerroelectric Ordered Water in Kaolinite Trenches using Monte CarloSimulations. J. Phys. Chem. A 2010. 114, 8396–8405. → pages 14, 105[91] Tunega, D.; Haberhauer, G.; Gerzabek, M. H.; Lischka, H. TheoreticalStudy of Adsorption Sites on the (001) Surfaces of 1:1 Clay Minerals.Langmuir 2002. 18, 139–147. → pages 14[92] Tunega, D.; Benco, L.; Haberhauer, G.; Gerzabek, M. H.; Lischka, H. AbInitio Molecular Dynamics Study of Adsorption Sites on the (001)Surfaces of 1:1 Dioctahedral Clay Minerals. J. Phys. Chem. B 2002. 106,11515–11525. → pages[93] Hu, X. L.; Michaelides, A. The Kaolinite (001) Polar Basal Plane. Surf.Sci. 2010. 604, 111–117. → pages 107, 108, 112[94] Hu, X. L.; Michaelides, A. Water on the Hydroxylated (001) Surface ofKaolinite: From Monomer Adsorption to a Flat 2D Wetting Layer. Surf.Sci. 2008. 602, 960–974. → pages 14[95] Hu, X. L.; Michaelides, A. Ice Formation on Kaolinite: Lattice Match orAmphoterism? Surf. Sci. 2007. 601, 5378–5381. → pages 107, 108, 112[96] Benco, L.; Tunega, D.; Hafner, J.; Lischka, H. Orientation of OH Groupsin Kaolinite and Dickite: Ab Initio Molecular Dynamics Study. Am.Mineral 2001. 86, 1057–1065. → pages 14[97] Pronk, S.; Pa´ll, S.; Schulz, R.; Larsson, P.; Bjelkmar, P.; Apostolov, R.;Shirts, M. R.; Smith, J. C.; Kasson, P. M.; van der Spool et al., D.GROMACS 4.5: A High-Throughput and Highly Parallel Open SourceMolecular Simulation Toolkit. Bioinformatics 2013. 29, 845–854. →pages 16[98] Frenkel, D.; Smit, B. Understanding Molecular Simulations. AcademicPress, San Diego, 2002. → pages 17, 18, 20, 21, 24, 25, 26157[99] van der Spoel, D.; Hess, L. B.; the GROMACS development team.GROMACS User Manual version 4.6.7. www.gromacs.org, 2014. →pages 18, 27, 28, 29[100] Abascal, J. L. F.; Sanz, E.; Ferna´ndez, R. G.; Vega, C. A Potential Modelfor the Study of Ices and Amorphous Water: TIP4P/Ice. J Chem Phys2005. 122, 234511. → pages 19, 29, 30, 40, 42, 85[101] Essman, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen,L. G. A Smooth Particle Mesh Ewald Method. J. Chem. Phys. 1995. 103,8577–8593. → pages 20, 22, 24[102] Darden, T.; York, D.; Pedersen, L. Particle Mesh Ewald: An Nlog(N)Method of Ewald Sums in Large Systems. J. Chem. Phys. 1993. 98,10089–10092. → pages[103] Lee, S. An FPGA Implementation of the Smooth Particle Mesh EwaldReciprocal Sum Compute Engine (RSCE). Master’s thesis, University ofToronto, 2005. → pages 22, 24[104] Nose´, S. A Unified Formulation of the Constant Temperature MolecularDynamics Methods. J. Chem. Phys. 1984. 81, 511–519. → pages 26[105] Hoover, W. G. Canonical Dynamics: Equilibrium Phase-SpaceDistributions. Phys. Rev. A. 1985. 31, 1695–1697. → pages 26[106] Bussi, G.; Donadio, D.; Parrinello, M. Canonical Sampling ThroughVelocity Rescaling. J. Chem. Phys. 2007. 126, 014101. → pages 27[107] Parrinello, M.; Rahman, A. Polymorphic Transitions in Single Crystals: ANew Molecular Dynamics Method. J. Apple. Phys. 1981. 52, 7182–7190.→ pages 27, 66[108] Nose´, S.; Klein, M. L. Constant Pressure Molecular Dynamics forMolecular Systems. Mol. Phys. 1983. 50, 1055–1076. → pages 27, 66[109] Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. G. E. M. LINCS: ALinear Constraint Slover for Molecular Simulations. J. Comput. Chem.1997. 18, 1463–1472. → pages 28158[110] Nada, H.; van der Earden, J. P. J. M. An Intermolecular Potential Modelfor the Simulation of Ice and Water Near the Melting Point: a six-siteModel of H2O. J Chem Phys 2003. 118, 7401–7413. → pages 29, 30, 40,42, 107[111] Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein,M. L. Comparison of Simple Potential Functions for Simulating LiquidWater. J. Chem. Phys. 1983. 79, 926–935. → pages 29, 85[112] Ferna´ndez, R. G.; Abascal, J. L. F.; Vega, C. The Melting Point of Ice Ihfor Common Water Models Calculated from Direct Coexistence of theSolid-Liquid Interface. J Chem Phys 2006. 124, 144506. → pages 29, 87[113] Abascal, J. L. F.; Ferna´ndez, R. G.; Vega, C.; Carignano, M. A. TheMelting Temperature of the Six Site Potential Model of Water. J ChemPhys 2006. 125, 166101. → pages 29, 87[114] S Grazˇulis, S.; Dasˇkevicˇ, A.; Merkys, A.; Chateigner, D.; Lutterotti, L.;Quiro´s, M.; Serebryanaya, N. R.; Moeck, P.; Downs, R. T.; et al., A. L. B.Crystallography Open Database (COD) - An Open-Access Collection ofCrystal Structures and Platform for World-Wide Collaboration. Nucleic.Acids. Res. 2011. 40, D420–D427. → pages 30[115] Downs, R. T.; Hall-Wallace, M. The American Mineralogist CrystalStructure Database. Am. Mineral. 2003. 88, 247–250. → pages[116] Berman, H. M.; Westbrook, J.; Feng, Z.; Gilliland, G.; Bhat, T. N.;Weissig, H.; Shindyalov, I. N.; Bourne, P. E. The Protein Data Bank.Nucleic. Acids. Res. 2000. 28, 235–242. → pages 30, 137[117] Hull, S.; Keen, D. A. Pressure-Induced Phase Transitions it AgCl, AgBr,and AgI. Phys Rev B 1999. 59, 750–761. → pages 30[118] Bish, D. L. Rietveld Refinement of the Kaolinite Structure at 1.5 K. ClaysClay Miner. 1993. 41, 738–744. → pages 30, 107[119] Prince, E.; Donnay, G.; Martin, R. F. Neutron Diffraction Refinement ofan Ordered Orthoclase Structure. Am. Mineral. 1973. 58, 500–507. →pages 30159[120] Saalfeld, H.; Wedde, M. Refinement of the Crystal Structure of Gibbsite,Al(OH)3. Z. Kristallogr. 1974. 139, 129–135. → pages 30, 127[121] Middleton, A. J.; Marshall, C. B.; Faucher, F.; Bar-Dolev, M.; Braslavsky,I.; Campbell, R. L.; Walker, V. K.; Davies, P. L. Antifreeze Protein fromFreeze-Tolerant Grass has a Beta-Roll with an Irregularly StructuredIce-Binding Site. J. Mol. Biol. 2012. 416, 713–724. → pages 30, 137, 146[122] Cygan, R. T.; Liang, J.-J.; Kalinichev, A. G. Molecular Models ofHydroxide, Oxyhydroxide, and Clay Phases and the Development of aGeneral Force Field. J. Phys. Chem. B 2004. 108, 1255–1266. → pages32, 112, 122[123] Zerze, G. H.; Mittal, J.; Best, R. B. Diffusive Dynamics of ContactFormation in Disordered Polypeptides. Phys. Rev. Lett. 2016. 116,068102. → pages 32[124] Duan, Y.; Wu, C.; Chowdhury, S.; Lee, M. C.; Xiong, G.; Zhang, W.;Yang, R.; Cieplak, P.; Luo, R.; Lee, T.; Caldwell, J.; Wang, J.; Kollman, P.A Point-Charge Force Field for Molecular Mechanics and Simulations ofProteins Based on Condensed-Phase Quantum Mechanical Calculations.J. Comput. Chem 2003. 24, 1999–2012. → pages 32, 139[125] Shelley, J. C.; Patey, G. N. Boundary Condition Effects in Simulations ofWater Confined Between Planar Walls. Mol. Phys. 1996. 88, 385–398. →pages 33[126] Hall, K.; Ashtari, M.; Cann, N. M. On Simulations of Complex Interfaces:Molecular Dynamics Simulations of Stationary Phases. J. Chem. Phys.2012. 136, 114705. → pages 33[127] Svishchev, I. M.; Kusalik, P. G. Electrofreezing of Liquid Water: AMicroscopic Perspective. J. Am. Chem. Soc. 1996. 188, 649–654. →pages 34[128] Humphrey, W.; Dalke, A.; Schulten, K. VMD: Visual MolecularDynamics. J. Mol. Graphics. 1996. 14, 33–38. → pages 35[129] Schro¨dinger, LLC. The PyMOL Molecular Graphics System, Version 1.8,2015. → pages 35160[130] Moore, E. B.; de la Llave, E.; Welke, K.; Scherlis, D. A.; Molinero, V.Freezing, Melting and Structure of Ice in a Hydrophilic Nanopore. Phys.Chem. Chem. Phys. 2010. 12, 4124–4134. → pages 35, 37, 45, 76[131] Ester, M.; Kriegel, H.-P.; Sander, J.; Xu, X. A Density-Based Algorithmfor Discovering Clusters in Large Spatial Databases with Noise. Proc. 2ndInt. Conf. on Knowledge Discovery and Data Mining 1996. 226–231. →pages 37, 76[132] Thangaraj, K.; Palanisamy, M.; Gobinathan, R.; Ramasamy, P.Dependence of Ice Nucleating Ability on Misfit. Journal of MaterialScience Letters 1986. 5, 326–328. → pages 44[133] Cox, S. J.; Kathmann, S. M.; Purton, J. A.; Gillan, M. J.; Michaelides, A.Non-Hexagonal Ice at Hexagonal Surfaces: The Role of LatticeMismatch. Phys. Chem. Chem. Phys 2012. 14, 7944–7949. → pages 44[134] Tasker, P. W. The Stability of Ionic Crystal Surfaces. J. Phys. C: SolidState Phys. 1979. 12, 4977–4984. → pages 62[135] Zielke, S. A.; Bertram, A. K.; Patey, G. N. A Molecular Mechanism of IceNucleation on Model AgI Surfaces. J. Phys. Chem. B 2015. DOI:10.1021/jp508601s. → pages 70, 172, 174[136] Mori, A.; Maruyama, M.; Furukawa, Y. Second-Order Expansion ofGibbs-Thomson Equation and Melting Point Depression of Ice Crystallite.J. Phys. Soc. Jpn. 1996. 65, 2742–2744. → pages 75[137] Pummer, B. G.; Budke, C.; Augustin-Bauditz, S.; Niedermeier, D.;Felgitsch, L.; Kampf, C. J.; Huber, R. G.; Liedl, K.; Loerting, T.;Moschen, T.; et al. Ice Nucleation by Water-Soluble Macromolecules.Atoms. Chem. Phys. 2015. 15, 4077–4091. → pages 76[138] Seo, M.; Jang, E.; Kim, K.; Choi, S.; Kim, J. S. UnderstandingAnisotropic Growth Behavior of Hexagonal Ice on a Molecular Scale: AMolecular Dynamics Simulation Study. J. Chem. Phys. 2012. 137,154503–1–8. → pages 81[139] Liu, J.; Nicholson, C. E.; Cooper, S. J. Direct Measurment of CriticalNucleus Size in Confined Volumes. Langmuir 2007. 23, 7286–7292. →pages 86161[140] Kra¨mer, B.; Hu¨bner, O.; Vortisch, H.; Wo¨ste, L.; Leisner, T.; Schwell, M.;Ru¨hl, E.; Baumga¨rtel, H. Homogeneous Nucleation Rates of SupercooledWater Measured in Single Levitated Microdroplets. J. Chem. Phys. 1999.111, 6521–6527. → pages 86[141] Malaspina, D. C.; di Lorenzo, A. J. B.; Pereyra, R. G.; Szleifer, I.;Carignano, M. A. The Water Supercooled Regime as Described by FourCommon Water Models. J. Chem. Phys. 2013. 139, 024506–1–8. →pages 87[142] Tenney, C. M.; Cygan, R. T. Molecular Simulation of Carbon Dioxide,Brine, and Clay Mineral Interactions and Determination of ContactAngles. Environ. Sci. Technol. 2014. 48, 2035–2042. → pages 94[143] Fenter, P.; Teng, H.; Geissbuhler, P.; Hanchar, J. M.; Nagy, K. L.; Sturchio,N. C. Atomic-Scale Structure of the Orthoclase (001)Water InterfaceMeasured with High-Resolution X-Ray Reflectivity. Geochim.Cosmochim. Ac. 2000. 64, 3663–3673. → pages 122, 127[144] Fenter, P.; Cheng, L.; Park, C.; Zhang, Z.; Sturchio, N. C. Structure of theOrthoclase (001)- and (010)-Water Interfaces by High-Resolution X-RayReflectivity. Geochim. Cosmochim. Ac. 2003. 67, 4267–4275. → pages122[145] Pedevilla, P.; Cox, S. J.; Slater, B.; Michaelides, A. Can Ice-LikeStructures Form on Non-Ice-Like Substrates? The Example of theKFeldspar Microcline. J. Phys. Chem. C 2016. 120, 6704–6713. → pages125[146] Zolles, T.; Burkart, J.; Hausler, T.; Pummer, B.; Hitzenberger, R.; Grothe,H. Identification of Ice Nucleation Active Sites on Feldspar Dust Particles.J. Phys. Chem. A 2015. 119, 2692–2700. → pages 126, 127[147] Nezbeda, I.; Moucka, F.; Smith, W. R. Recent Progress in MolecularSimulation of Aqueous Electrolytes: Force Fields, Chemical Potentialsand Solubility. Mol. Phys. 2016. 114, 1665–1690. → pages 126[148] Kulkami, G.; Sanders, C.; Zhang, K.; Liu, X.; c Zhao. Ice Nucleation ofBare and Sulfuric Acid-Coated Mineral Dust Particles and Implication forCloud Properties. J. Geophys. Res. Atmos. 2014. 119, 9993–10011. →pages 127162[149] Crooks, G. E.; Hon, G.; Chandonia, J. M.; Brenner, S. E. WebLogo: ASequence Logo Generator. Genome Res 2004. 14, 1188–1190. → pages135[150] Graether, S. P.; Jia, Z. Modeling Pseudomonas syringae Ice-NucleationProtein as a β -Helical Protein. Biophys. J. 2001. 80, 1169–1173. → pages135[151] Pandey, R.; Usui, K.; Livingstone, R. A.; Fischer, S. A.; Pfaendtner, J.;Backus, E. H. G.; Nagata, Y.; Fro¨hlich-Nowoisky, J.; Schmu¨ser, L.; Mauri,S.; et al. Ice-Nucleating Bacteria Control the Order and Dynamics ofInterfacial Water. Sci. Adv. 2016. 2, e1501630. → pages 135[152] Vrielink, A. S. O.; Aloi, A.; Olijve, L. L. C.; Voets, I. K. Interaction of IceBinding Proteins with Ice, Water and Ions. Biointerphases 2016. 11,018906. → pages 137[153] Middleton, A. J.; Brown, A. M.; Davies, P. L.; Walker, V. K. Identificationof the Ice-Binding Face of a Plant Antifreeze Protein. FEBS Lett 2009.583, 815–819. → pages 137163Appendix AAppendix to Chapter 3A.1 Results for the TIP4P/Ice modelIn general TIP4P/Ice behaves much like the six-site model. Overall, the densityprofiles for this model (Figures. A.2, A.3, and A.4) are similar to the profiles re-ported for six-site water. One notable difference is the appearance of an additionalsurface peak in the liquid phase profiles (Figures. A.2B, A.3B, and A.4B), and toa much lesser extent in those of ice (Figures. A.2A, A.3A, and A.4A). Evidently,liquid TIP4P/Ice has three preferred distances from the surface that are reduced totwo when the model freezes and enforces an ice structure.164ice molecules01,0002,0003,0004,000time at 250K (ns)0 20 40 60 80 1001Figure A.1: The number of ice molecules as a function of time for theTIP4P/Ice water model detected by the CHILL algorithm in simu-lations on three AgI surfaces: Ag exposed β -AgI(0001) (red), Agexposed γ-AgI(001)( blue), and Ag exposed γ-AgI(111) (green).The solid lines and dotted lines represent the number of Ic and Ihmolecules, respectively.Table A.1: List of simulations performed using the six-site model of waterat 270 K (supercooled by 20 K)AgI Phase Surface Run time (ns) Ice obtained?βAg exposed (0001) 210 YesI exposed (0001) 508 NoAg+I exposed (101¯0) 313 NoγAg exposed (111) 224 YesI exposed (111) 209 NoAg exposed (001) 98 YesI exposed (001) 198 No165 ρ/ρ o02468101214distance from surface (nm)0 0.5 1ABC1Figure A.2: Density profiles of oxygen (red) and hydrogen (black) atoms forthe TIP4P/Ice model near β -AgI(0001) surfaces. Panels A and B arefor ice (supercooled by 20 K) and liquid water (5 K above the meltingpoint), respectively, near the silver exposed face. Panel C is for liquidwater (supercooled by 20 K) near the iodide exposed face. On thedistance scale, zero (vertical green line) coincides with the positionof the centre of the outermost ion of the surface. ρo is the averagedensity far from the surface. Note the shifted scale in Panels A and B.Table A.2: List of simulations performed using the TIP4P/Ice water modelat 250 K (supercooled by 20 K)AgI Phase Surface Run time (ns) Ice obtained?βAg exposed (0001) 209 YesI exposed (0001) 509 NoAg+I exposed (101¯0) 313 NoγAg exposed (111) 209 YesI exposed (111) 209 NoAg exposed (001) 98 YesI exposed (001) 198 No166 ρ/ρ o02468101214distance from surface (nm)0 0.2 0.4 0.6ABC1Figure A.3: Density profiles of oxygen (red) and hydrogen (black) atoms forthe TIP4P/Ice model near γ-AgI(001) surfaces. Panels A and B are forice (supercooled by 20 K) and liquid water (10 K above the meltingpoint), respectively, near the silver exposed face. Panel C is for liquidwater (supercooled by 20 K) near the iodide exposed face. On thedistance scale, zero (vertical green line) coincides with the positionof the centre of the outermost ion of the surface. ρo is the averagedensity far from the surface. Note the shifted scale in Panels A and B.Table A.3: List of simulations performed to determine the highest tempera-ture at which ice will nucleate on the Ag exposed β -AgI(0001) faceWater model Temperature (K) Run time (ns) Ice obtained? Notessix-site280 49 Yes281 155 Yes froze after 76 ns282 255 No283 55 No285 55 NoTIP4P/Ice260 209 Yes271 249 Yes froze after 149 ns272 249 No273 49 No275 209 No167 ρ/ρ o02468101214distance from surface (nm)0 0.5 1ABC1Figure A.4: Density profiles of oxygen (red) and hydrogen (black) atoms forthe TIP4P/Ice model near γ-AgI(111) surfaces. Panels A and B arefor ice (supercooled by 20 K) and liquid water (5 K above the meltingpoint), respectively, near the silver exposed face. Panel C is for liquidwater (supercooled by 20 K) near the iodide exposed face. On thedistance scale, zero (vertical green line) coincides with the positionof the centre of the outermost ion of the surface. ρo is the averagedensity far from the surface. Note the shifted scale in Panels A and B.168 ρ/ρ o00.511.522.53distance from surface (nm)0 0.5 1 1.51Figure A.5: Density profiles of oxygen (red) and hydrogen (black) atomsfor the TIP4P/Ice model near the Ag+I exposed β -AgI(101¯0) surface(supercooled by 20 K). On the distance scale, zero (vertical green line)coincides with the position of the centre of the outermost ion of thesurface. ρo is the average density far from the surface.Table A.4: List of simulations performed using non-standard charges. Thesix-site simulation with charges of ±0.8 e on the β -AgI(0001) surfaceonly froze after 473 ns.Water model Ion exposed Charge (e) Ice obtained?/Run time (ns)β -AgI(0001) γ-AgI(001)six-siteI 0.0 No/255 No/317Ag 0.0 Yes/55 No/317Ag 0.2 Yes/55 Yes/27Ag 0.4 Yes/55 Yes/27Ag 0.8 Yes/555 No/27Ag 1.0 No/655 —TIP4P/IceI 0.0 No/98 No/98Ag 0.0 Yes/98 No/98Ag 0.2 Yes/98 Yes/27Ag 0.4 Yes/98 Yes/27Ag 0.8 No/298 No/27Ag 1.0 No/198 —169 ρ/ρ o02468distance from surface (nm)0 0.5 1AB1Figure A.6: Density profiles of oxygen (red) and hydrogen (black) atomsfor the TIP4P/Ice model near the Ag exposed β -AgI(0001) surface .Panel A is for neutral AgI (5 K above the melting point) and PanelB is for AgI with charges of ± 1.0e (supercooled by 20 K). On thedistance scale, zero (vertical green line) coincides with the positionof the centre of the outermost ion of the surface. ρo is the averagedensity far from the surface.170 ρ/ρ o01234567distance from surface (nm)0 0.5 1AB1Figure A.7: Density profiles of oxygen (red) and hydrogen (black) atoms forthe TIP4P/Ice model near the Ag exposed γ-AgI(001) surface (super-cooled by 20 K). Panel A is for neutral AgI and Panel B is for AgIwith charges of ± 1.0e. On the distance scale, zero (vertical greenline) coincides with the position of the centre of the outermost ion ofthe surface. ρo is the average density far from the surface.171Appendix BAppendix to Chapter 4The simulations used in this study are listed in Table B.1 for the six-site model,and in Table B.3 for the TIP4P/2005 model. For each simulation, the system,temperature, and length of the simulation are given. Whether or not ice nucleatedis indicated in the last column.B.1 Comparison of force fieldsIn addition to the Lennard-Jones and Coulombic terms used in our present andprevious[135] simulations, the original force field proposed by Hale and Kiefer[35]for a water molecule interacting with an AgI surface included a polarization inter-action of the general form172uind =−0.5∑n∑mαwqnqm(ro− rn)(ro− rm)|ro− rn|3|ro− rm|3 −0.5∑m ∑i ∑jαmqiq j(ro− r i)(ro− r j)|ro− rm|3|ro− rm|3 ,(B.1)where ε0 is the permittivity of free space, αw and αm are the polarizabilities ofwater and of Ag or I ions, respectively, Qk represents the partial charge at a sitek on the AgI lattice or on the water molecule, rk is the position vector of any sitek, and ro is the position vector of the oxygen atom of the water molecule. Thesums on n and m are over all Ag and I ions in the AgI lattice, and the sums oni and j are over all charge sites in the water model. The first term in Eq. B.1represents the interaction of the AgI lattice with the induced dipole moment of thewater molecule, and the second term is the interaction of the water molecule withthe dipole moment induced in a Ag or I ion.In their calculations Hale and Kiefer[35] actually used an approximate form ofthe AgI-water induction potential, retaining only the m = n terms in the first termof Eq. B.1. They note that this is a good approximation because the m 6= n termscontribute only about 5% to the total potential. This simplification means thatthe water-molecule-AgI-surface interaction can be written as a sum of pairwiseinteractionsuIND =∑mu˜IND, (B.2)where the sum on m is over all Ag and I ions, and u˜IND depends on the water-iondistance and on the orientation of the water molecule. Note that the induction173potential is still not pairwise additive with respect to the multiple charge sitesin the water model [note the sums on i and j in the second term of Eq. B.1],which means that u˜IND cannot be efficiently included in parallel MD algorithms,however, it does allow us to meaningfully compare total Ag-water and I-waterpair potentials both with and without induction effects.Results obtained for the six-site water model for three orientations of the wa-ter molecule (see Fig. B.1) are shown in Fig. B.2. We note that results for theTIP4P/ice model (not shown) are very similar. In Fig. B.2 four curves are shownfor each orientation of the water molecule. The Hale-Kiefer (HK) potentials withAg and I partial charges of ±0.6e and including u˜IND are represented by greencurves. HK potentials neglecting u˜IND, as used in the present paper, are shownas black curves. Additionally, we include curves for potentials that neglect u˜INDbut have Ag and I partial charges of zero (blue curves) and ±0.8e (red curves).In an earlier paper [135] we showed that AgI slabs modelled with partial chargesranging from zero to ±0.8e all nucleated ice in a similar manner. From Fig. B.2,we see that while the HK potential does have a significant contribution from u˜IND(particularly for attractive orientations) it does not fall outside the range of poten-tials for which we have observed ice nucleation. In particular, for the importantorientations of the water molecule where the interactions are strongly attractive,the HK potential is comparable with the potential obtained with partial chargesof ±0.8e. Therefore, based on our observation of ice nucleation for a wide rangeof Coulombic interactions, including the three cases shown in Fig. B.2, we donot believe that using the full HK potential (if that were computationally feasible)174would qualitatively alter our results or conclusions.175Table B.1: A list of simulations performed with the six-site model. Systemsmarked with a ‘P’ were conducted at 1 bar pressure and a ‘0’ indicatesthe Ag and I atoms were neutral.System ∆Ts (K) Time (ns) Nucleated ice?A 31 50 NoA 29 50 NoA 34 150 NoA 39 102 YesA 36 127 YesA 37 100 YesA 35 150 YesA 35 150 YesA 34 150 NoAP 35 300 YesB 29 98 YesB 24 98 NoB 26 150 YesB 27 98 YesB 25 39 NoB 26 248 YesB 25 39 NoB 21 150 NoBP 26 175 YesC 19 100 NoC 24 48 YesC 21 98 NoC 22 78 YesC 22 78 YesC 21 98 NoC 22 100 YesC 17 100 NoC 14 100 NoCP 22 125 YesC0 22 100 YesC0 21 100 YesC2 29 48 YesC2 29 48 Yes176Table B.2: Continuation of Table B.1.System ∆Ts (K) Time (ns) Nucleated ice?C2 24 48 YesC2 22 48 YesC2 19 48 NoC2 20 98 NoC2 22 100 YesC2 21 100 YesC2 20 100 NoC2 19 100 NoC2 21 100 YesD 14 98 NoD 16 148 YesD 15 150 NoD 15 150 NoD 16 123 YesD 24 50 YesD 29 50 YesD 24 50 YesD 29 50 YesD 13 100 NoD 9 100 NoDP 18 75 YesD0 16 150 YesD0 18 75 YesDr 16 115 YesDr 15 150 NoDr 16 113 YesDr 16 183 YesDr 15 150 NoE 13 150 NoE 12 150 NoE 14 175 YesE 15 100 YesE 9 150 NoE 14 150 YesE 13 150 No177Table B.3: A list of simulations performed with the TIP4P/2005 modelSystem ∆Ts (K) Time (ns) Nucleated ice?D 22.5 300 YesD 28.5 300 YesD 25.5 300 YesD 30.5 200 YesD 20.5 600 YesD 21.5 500 YesD 19.5 500 YesD 18.5 500 NoD 18.5 600 YesD 19.5 500 NoE 20.5 400 YesE 15.5 300 NoE 25.5 300 YesE 17.5 600 NoE 18.5 600 YesE 19.5 100 YesE 17.5 600 YesE 18.5 500 Yes178ABCrFigure B.1: Orientations of water models tested to produce interaction ener-gies. Purple circles represent the ion, either Ag or I, and red and blackcircles represent the oxygen and hydrogen atoms, respectively, of wa-ter. Orientation A has the dipole moment of water pointing towardthe ion. In B, the dipole moment of water is pointing away from theion, and in C a hydrogen bond is pointing toward the ion. The arrowlabeled ‘r’ shows the direction the water was moved in the calculationof the interaction energies.179Figure B.2: Interaction energies between an Ag or I ion and a six-site watermodel in the three orientations shown is Figure B.1. The title of eachplot lists the orientation and the ion involved. ZBP curves are from thepotential used in this study. ZBP 0.8e and ZBP 0.0e are the presentpotential using charges of± 0.8e and 0.0 on the AgI ions, respectively.HK is the potential[35] from which the ZBP potential is derived fromand includes polarization terms.180Appendix CAppendix to Chapter 5Table C.1: Rigid-kaolinite simulations.Model Time (ns) Temperature (K) Ice Nucleated?TIP4P/Ice 100 240 NoTIP4P/Ice 500 230 NoTIP4P/Ice 500 230 Nosix-site 50 240 Nosix-site 500 240 Nosix-site 500 240 No181Table C.2: Free-H simulations.Model Time (ns) Temperature (K) Ice Nucleated?TIP4P/Ice 198 240 NoTIP4P/Ice 398 230 YesTIP4P/Ice 398 230 YesTIP4P/Ice 327 230 YesTIP4P/Ice 377 230 Yessix-site 589 240 Yessix-site 300 230 Nosix-site 300 240 Nosix-site 700 240 YesTable C.3: Free-OH simulationsModel Time (ns) Temperature (K) Ice Nucleated?TIP4P/Ice 500 230 NoTIP4P/Ice 500 240 NoTIP4P/Ice 500 230 Nosix-site 500 240 Nosix-site 500 250 Nosix-site 500 240 No182Table C.4: Free-OH simulations started from free-H simulations, as de-scribed in Chapter 5. TIP4P/Ice simulations were at 230 K and six-sitesimulations at 240 K. Divisions in table separate simulations startedfrom different Free-H simulations.Model Time (ns) Ice Nucleated? Starting time (ns)TIP4P/Ice 100 No 85TIP4P/Ice 100 No 165TIP4P/Ice 200 Yes 185TIP4P/Ice 100 Yes 200TIP4P/Ice 75 Yes 250TIP4P/Ice 100 No 110TIP4P/Ice 200 Yes 190TIP4P/Ice 200 Yes 205six-site 100 No 150six-site 100 Yes 225six-site 100 No 270six-site 125 Yes 300six-site 100 Yes 340six-site 100 Yes 400183Table C.5: Ice-prepared surface simulations. Simulations marked with anasterisk used ice-prepared surfaces produced by the other model of wa-ter.Model Time (ns) Temperature (K) Ice Nucleated?TIP4P/Ice 171 240 YesTIP4P/Ice 50 260 NoTIP4P/Ice 150 255 YesTIP4P/Ice 171 240 YesTIP4P/Ice 50 250 YesTIP4P/Ice 150 255 YesTIP4P/Ice* 200 255 YesTIP4P/Ice* 200 255 Yessix-site 346 240 Yessix-site 50 280 Nosix-site 50 260 Yessix-site 138 250 Yessix-site 50 280 Nosix-site 50 260 Yessix-site* 50 260 Yessix-site* 260 260 YesTable C.6: Liquid-prepared surface simulations.Model Time (ns) Temperature (K) Ice Nucleated?TIP4P/Ice 150 240 NoTIP4P/Ice 700 230 NoTIP4P/Ice 150 240 NoTIP4P/Ice 500 230 Nosix-site 50 260 Nosix-site 500 240 Nosix-site 500 240 No184Table C.7: “Ideal” surface simulations.Model Time (ns) Temperature (K) Ice Nucleated?TIP4P/Ice 600 230 YesTIP4P/Ice 600 230 Yessix-site 650 240 Yessix-site 502 240 NoTable C.8: Vacuum-prepared surface simulations.Model Time (ns) Temperature (K) Ice Nucleated?TIP4P/Ice 300 230 NoTIP4P/Ice 500 230 NoTIP4P/Ice 500 230 Nosix-site 300 240 Nosix-site 500 240 Nosix-site 500 240 NoTable C.9: Si-surface simulations. The simulation marked with an asteriskis the large Si-surface simulation.Model Time (ns) Temperature (K) Ice Nucleated?TIP4P/Ice 100 240 NoTIP4P/Ice 1100 230 YesTIP4P/Ice 1100 230 YesTIP4P/Ice 700 220 NoTIP4P/Ice 200 250 Nosix-site 100 240 Nosix-site 700 240 Yessix-site 800 240 Yessix-site 200 250 Yessix-site* 373 250 Yes185

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0314173/manifest

Comment

Related Items