UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Determination of melting point trends of model salts by molecular dynamics simulations Lindenberg, Erin 2014

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

Item Metadata


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

Full Text

Determination of melting point trends of model saltsby molecular dynamics simulationsbyErin LindenbergA THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFDOCTOR OF PHILOSOPHYinThe Faculty of Graduate and Postdoctoral Studies(Chemistry)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)December 2014c© Erin Lindenberg 2014AbstractWe examine the melting point trends across sets of coarse grain model saltsusing NPT molecular dynamics simulations. The melting point trends are estab-lished relative to a charge-centered, size-symmetric salt that is closely akin to therestricted primitive model. Two of the common features of ionic liquids, namelysize asymmetry and a distributed cation charge, are systematically varied in a setof model salts. We find that redistributing the cation charge in salts with size-symmetric, monovalent, spherical ions can reduce the melting temperature by upto 50% compared to the charge-centered case. Displacing the charge from the ioncenter reduces the enthalpy of the liquid more than that of the solid resulting in alower melting point. We consider two sets of size-asymmetric salts with size ratiosup to 3:1 using different length scales; the melting point trends are different in eachset, but within each set we find salts that achieve a melting point reduction of over60% from the charge-centered, size-symmetric case. The lowest melting point rangewe find is between 450 K and 500 K. We find diversity in the solid phase structures.For all size ratios with small cation charge displacements, the salts crystallize withorientationally disordered cations. For equal-sized ions, once the cation charge ismoved far enough off-center, the salts become trapped in glassy states upon coolingand we find an underlying crystal structure (space group 111) that features orien-tationally ordered ion pairs. The salts with large size ratios and large cation chargedisplacements achieve the lowest melting points and also show premelting transi-iitions at lower temperatures (two as low as 300 K). We find two types of premeltingbehaviour; some salts exhibit a fast ion conductor phase, where the smaller anionsmove through a face-centered cubic (fcc) cation lattice, whereas other salts have aplastic crystal phase composed of ion pairs rotating on an fcc lattice.iiiPrefaceThe research presented in this thesis are or will be co-authored publicationsof E. K. Lindenberg and Dr. G. N. Patey. Parts of Chapters 1, 2, 3, and 6 have beenpublished as a journal article [E. K. Lindenberg and G. N. Patey, “How distributedcharge reduces the melting points of model ionic salts,” J. Chem. Phys. 140,104504 (2014)]. Parts of Chapters 1, 2, 4, and 6 are in preparation for submission asa journal article. The research and data analysis was performed by E. K. Lindenbergwith guidance and suggestions from Dr. G. N. Patey. The manuscripts were writtenby E. K. Lindenberg with revisions and additions by Dr. G. N. Patey.The work in this thesis builds on earlier projects done by Dr. H. V. Spohras a doctoral student and post-doctoral researcher in Dr. G. N. Patey’s group. Dr.H. V. Spohr’s doctoral work focused on the liquid phase dynamics of similar modelsalts. Part of Dr. H. V. Spohr’s post-doctoral research was a preliminary study ofthe melting point trends of a set of charge- and size-symmetric salts.The model salts considered in thesis include and expand upon the ion modelgeometries from Dr. H. V. Spohr’s doctoral work. Some of the ion model param-eters were set by Dr. H. V. Spohr and Dr. G. N. Patey. The methods used todetermine melting point trends in this thesis differ from those used in Dr. H. V.Spohr’s preliminary melting point work. The final research design, data collectionand analysis, and manuscript composition was done by E. K. Lindenberg with guid-ance and suggestions from Dr. G. N. Patey, and was carried out independently fromthe initial work of Dr. H. V. Spohr.ivTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiiList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixList of Symbols . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiv1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Ionic Liquids . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Computational Studies of ILs . . . . . . . . . . . . . . . . . . . . . . 61.3 The Melting Points of Ionic Liquids . . . . . . . . . . . . . . . . . . 71.4 Experimental Melting Points . . . . . . . . . . . . . . . . . . . . . . 101.5 Simulation Methods for Estimating Melting Temperatures . . . . . . 122 Models and Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142.1 Ion Topologies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142.2 Ion Interaction Parameters . . . . . . . . . . . . . . . . . . . . . . . 242.3 Ion Masses . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 272.4 Charge Arm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 282.5 Minimum Potential Energies of Isolated Ion Pairs, U IPiso . . . . . . . . 292.6 Outline of the Salts Studied in This Thesis . . . . . . . . . . . . . . 30v2.7 Hysteresis Method for Estimating Melting Temperatures . . . . . . . 312.8 Two-Phase Simulation Method for Estimating Melting Temperatures 342.9 Simulation Details . . . . . . . . . . . . . . . . . . . . . . . . . . . . 372.10 Crystal Structures . . . . . . . . . . . . . . . . . . . . . . . . . . . . 403 Cation Charge Distribution: Size Symmetric 2L - 1C Salts (41Salts) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 423.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 423.2 Melting Point Results . . . . . . . . . . . . . . . . . . . . . . . . . . 433.3 Melting Point Trends . . . . . . . . . . . . . . . . . . . . . . . . . . 593.4 2L - 1C Salt Properties . . . . . . . . . . . . . . . . . . . . . . . . . 643.5 Chapter Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . 784 Ion Size Ratios: 2L - 1C (105 Salts) . . . . . . . . . . . . . . . . . . 804.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 804.2 Size Set A (Constrained σ±) . . . . . . . . . . . . . . . . . . . . . . 814.3 Size Set B (Unconstrained σ±) . . . . . . . . . . . . . . . . . . . . . 954.4 Solid Structure and Ion Dynamics . . . . . . . . . . . . . . . . . . . 1034.5 Chapter Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1105 Cation Charge Geometry: Size Symmetric 3s - 1C Salts (15 Salts) 1115.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1115.2 Melting Point Results . . . . . . . . . . . . . . . . . . . . . . . . . . 1125.3 Chapter Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1186 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1196.1 Analysis Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1196.2 Connection with Experimental Work . . . . . . . . . . . . . . . . . . 1226.3 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130viAppendicesA Configurational Calculations . . . . . . . . . . . . . . . . . . . . . . . 142A.1 Simulation Cell . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143A.2 Particle Representations . . . . . . . . . . . . . . . . . . . . . . . . . 145A.3 Interaction Potentials, Forces, Pressure . . . . . . . . . . . . . . . . . 147B Molecular Dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155B.1 Time Evolution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155B.2 Simulation Constraints . . . . . . . . . . . . . . . . . . . . . . . . . 158B.2.1 Temperature — Particle Motion . . . . . . . . . . . . . . . . 158B.2.2 Pressure — Simulation Cell . . . . . . . . . . . . . . . . . . . 160C Genetic Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 162C.1 Initialization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 163C.2 Configuration Checks . . . . . . . . . . . . . . . . . . . . . . . . . . 164C.3 Configuration Evolution . . . . . . . . . . . . . . . . . . . . . . . . . 165C.3.1 Diversity Operators . . . . . . . . . . . . . . . . . . . . . . . 166C.3.2 Refinement Operators . . . . . . . . . . . . . . . . . . . . . . 169C.4 Termination . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 170C.5 GA Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 171D Temperature Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . 172E Crystal Structures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 175E.1 111n Structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 175E.2 Var-CsCl Structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . 177E.3 Var-111n Structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . 178F Input Files . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 179viiList of Tables2.1 Parameters of the Size-Asymmetric Salts . . . . . . . . . . . . . . . . 212.2 Rcut Values for the Size-Asymmetric Salts . . . . . . . . . . . . . . . 232.3 The Mass Distribution for the 2L and 3s Ions . . . . . . . . . . . . . 272.4 Salts Studied in Each Chapter of This Thesis . . . . . . . . . . . . . 313.1 Properties of the Spontaneously Crystallizing 2L - 1C Salts . . . . . . 493.2 Properties of the “Glassy” 2L - 1C Salts . . . . . . . . . . . . . . . . 593.3 Phases of 2L67 - 1C Salt Subset at 500 K and 1100 K . . . . . . . . . 694.1 Properties of the Set A Salts . . . . . . . . . . . . . . . . . . . . . . . 904.2 Properties of the Set B Salts . . . . . . . . . . . . . . . . . . . . . . . 965.1 Properties of the 3s - 1C Salts . . . . . . . . . . . . . . . . . . . . . . 1156.1 Melting Points of Alkylammonium Bromide Salts . . . . . . . . . . . 123B.1 Predictor-Corrector Numerical Coefficients . . . . . . . . . . . . . . . 158D.1 Melting Temperatures of the Crystallizing 2L - 1C Salts . . . . . . . . 172D.2 Melting Temperatures of the “Glassy” 2L - 1C Salts . . . . . . . . . . 174E.1 Unit Cell of the 111n Crystal Structure . . . . . . . . . . . . . . . . . 176viiiList of Figures1.1 Common Ions Found in Ionic Liquids . . . . . . . . . . . . . . . . . . 22.1 Diagram of the 2L Model Ion . . . . . . . . . . . . . . . . . . . . . . 162.2 Diagram of the 3s(120◦) Model Cation . . . . . . . . . . . . . . . . . 172.3 Diagrams of the Set of 3s Ions . . . . . . . . . . . . . . . . . . . . . . 192.4 Diagram of the Size-Asymmetric Salts . . . . . . . . . . . . . . . . . 202.5 3D Potential Surface Plots for Select 2L67 - 1C Salts . . . . . . . . . 262.6 Two-Phase Simulation Initial Configuration Snapshots . . . . . . . . 363.1 Hysteresis Plot of the 1C - 1C Salt . . . . . . . . . . . . . . . . . . . 443.2 Hysteresis Plots of a Subset of 2L - 1C Salts . . . . . . . . . . . . . . 473.3 Illustration of the 111n Crystal Structure . . . . . . . . . . . . . . . . 543.4 Hysteresis Plot of the 2L67-20 - 1C Salt . . . . . . . . . . . . . . . . 553.5 Hysteresis Plot of the 2L100-16 - 1C . . . . . . . . . . . . . . . . . . 563.6 Configurational Snapshots of the 2L100-16 - 1C Salt at T = 700 K . . 573.7 Melting Temperature Summary for the A100 2L - 1C Salts . . . . . . 603.8 Impact of Each 2L Ion Parameter on Melting Point . . . . . . . . . . 613.9 Plot of Scaled Melting Points Against Normalized Charge Arm . . . . 623.10 Average Enthalpies of 2L - 1C Salts at T = 1100 K . . . . . . . . . . 663.11 Average Enthalpies of 2L67 and 2L100 - 1C Salts at T = 600 K . . . 683.12 Subset of Radial Distribution Functions at T = 500 K . . . . . . . . . 713.13 Subset of Radial Distribution Functions at T = 1100 K . . . . . . . . 723.14 Subset of Velocity Autocorrelation Functions at T = 1100 K . . . . . 743.15 Subset of Rotational Correlation Functions at T = 500 K and T =1100 K . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 764.1 Hysteresis Plot of the A300 1C - 1C Salt . . . . . . . . . . . . . . . . 83ix4.2 Melting Temperature Summary for the 1C - 1C Salts in Set A . . . . 854.3 Hysteresis Plots of Four 2L67 - 1C Salts in Set A . . . . . . . . . . . 874.4 Melting Temperature Summary for the Set A Salts . . . . . . . . . . 934.5 Impact of d on Melting Points of Set A Salts . . . . . . . . . . . . . . 944.6 Melting Temperature Summary of the Set B Salts . . . . . . . . . . . 994.7 Plot of Scaled Melting Points Against Normalized Charge Arm . . . . 1024.8 Mean Square Displacements for the A300 Salts at T = 550 K . . . . . 1044.9 Anion Trajectories in Three A300 Salts at T = 550 K . . . . . . . . . 1064.10 Radial Distribution Functions of the A300 Salts at T = 550 K . . . . 1095.1 Hysteresis Plots of a Subset of 3s - 1C Salts . . . . . . . . . . . . . . 1135.2 Melting Temperature Summary of the 3s - 1C Salts . . . . . . . . . . 1176.1 Plot of the Normal Melting Points of Substituted Ammonium Bro-mide Salts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125A.1 Structures of MD and GA Programs . . . . . . . . . . . . . . . . . . 144E.1 Findsym Program Output on the 111n Unit Cell . . . . . . . . . . . . 175E.2 Var-CsCl Structure of the A100 2L67-22 - 1C Salt . . . . . . . . . . . 177E.3 Var-111n Structure for the B133 2L67-22 - 1C Salt . . . . . . . . . . . 178F.1 Example .pbs File . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 181F.2 Example .gro File . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 184F.3 Example .top File . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 185F.4 Example .mdp File . . . . . . . . . . . . . . . . . . . . . . . . . . . . 186xList of Symbolsµi(t) Unit Vector From CM to CC on 2L Cation i at Time tωli Angular Velocity of Particle i in the Local Particle Frameτ li Torque on Particle i in the Local Particle FrameA Rotation Matrix to Transform Coordinates Between Frames of Referencefi Force on Particle ig Metric Tensor Defining the Simulation Cellh Matrix that Defines the Conventional Simulation CellQ Quaternion Representation of a Particle’s OrientationrFi Fractional Coordinates of Particle i’s Center Within the Simulation Cellri Cartesian Coordinates of the Center (LJ Site) of Particle iriα Cartesian Coordinates of Interaction Site α on Particle isα Vector From the Particle Center to Interaction Site α in the ParticleCoordinate Frame Permittivity of the Mediumκ Ewald Summation Parameterρ∗ Reduced Number Density, Nσ3/Vσ The LJ Length Parameter (Particle Diameter)σ+ The LJ Diameter of the Cationσ− The LJ Diameter of the Anionσ± The LJ Distance From the Cation Center to the Anion Centerθ Angle Between Vectors Drawn From the Ion Center to the Two Off-Center Sites in a 3s Ionxiεij Lennard-Jones Energy Parameterξ Friction Coefficient that Adjusts Particle Motions to Maintain the Tem-peratureD Diffusion Coefficientd The Distance Between the Ion Center (LJ Site) and Off-center ChargeSite in Units of 0.01 nmLc Normalized Charge Armm MassN Number of ParticlesP PressureQ Parameter Determining the Thermostat Responsivity (Scalar)qion Total Ion ChargeRion Ion RadiusRcut Cutoff Distance for Computing Short-range Interactionsriαjβ The Distance Between Site α on Ion i and Site β on Ion jrij The Distance Between LJ Sites on Ions i and jT+ Highest Temperature of the Superheated Solid Before MeltingT− Lowest Temperature of the Supercooled Liquid Before FreezingTg Glass Transition Temperature Defined as the Temperature Where theDiffusion Coefficient Falls Below 10−11 m2/s on the Cooling CurveTl Lowest Temperature that Evolved into a Liquid in the Two-phase Sim-ulationsTm Melting Temperature Determined from the Hysteresis MethodTs Highest Temperature that Evolved into a Solid in the Two-phase Simu-lationsuij Pair Potential Between Particles i and jV Volume of the Simulation CellxiiW Parameter Determining the Barostat Responsivity (Scalar)1C A Spherical Ion with One Interaction Site at its Center2L An Ion with Two Interaction Sites in a Linear Arrangement2L100 A Two-site Ion with all of the Unit Charge on the Off-center InteractionSite2L33 A Two-site Ion with 1/3 of the Unit Charge on the Off-center InteractionSite2L50 A Two-site Ion with 1/2 of the Unit Charge on the Off-center InteractionSite2L67 A Two-site Ion with 2/3 of the Unit Charge on the Off-center InteractionSite3s Represents a Set of Ions, Each with Three Interaction SitesAC Anion Off-center Charge SiteAM Anion Center (LJ Site)CC Cation Off-center Charge SiteCM Cation Center (LJ Site)CsCl Two Interpenetrating Simple Cubic Lattices of Cations and AnionsIon Center Location of the Lennard-Jones Interaction Sitemsd Mean Squared DisplacementMSDS Material Safety Data SheetNaCl Two Interpenetrating Face Centered Cubic Lattices of Cations and An-ionsOD-CsCl Orientationally Disordered Cesium Chloride Crystal StructureRPM Restricted Primitive Modelxiiifor my familyxivChapter 1Introduction ∗“Because ionic liquids have complex molecular structure, the corre-lation of the properties displayed by these materials to their structure maybe very complicated; there is then a need to rationalize the behaviour of thisclass of liquids in a general picture, which is able to keep track of theiressential characteristics.”– Marco Malvaldi and Cinzia ChiappeJ. Phys.: Condens. Matter 20 035108 (2008)1.1 Ionic LiquidsIonic liquids (ILs) are salts usually composed of molecular ions that existin the liquid phase under ambient, or near ambient, conditions. By convention,ILs are defined as salts with melting points or glass transition temperatures be-low 373 K,1 which differentiates them from conventional molten salts that melt atmuch higher temperatures.2 Sodium chloride and cesium chloride, for example, havenormal melting points of 1074 K and 918 K, respectively.Ammonium-based ionic liquids have been known for over a century, with 2-aminoethanol nitrate reported in 1888 by S. Gabriel and J. Weiner, and ethylaminenitrate reported in 1914 by P. Walden.3,4 It was Walden who provided the initialdistinction of ILs as water-free salts which melt at relatively low temperatures, up∗Parts of this chapter have been published. E. K. Lindenberg and G. N. Patey, “How distributedcharge reduces the melting points of model ionic salts,” J. Chem. Phys. 140, 104504 (2014).1to about 373 K, which has become the accepted definition.3,4It is estimated that there are up to 106 possible binary ionic liquids, whichincreases to about 1018 if ternary systems are included.3,5 Over the last 50 years,research on ionic liquids has primarily focused on cation variations; tetraalkylammo-nium cations and organic heterocyclic cations such as imidazolium, pyridinium, andpyrrolidinium structures are common. Anions typically tend toward the inorganicside, such as halides, nitrate, tetrafluoroborate, and hexafluorophosphate, but arebecoming more exotic. Some of the common IL ions are shown in Fig. 1.1.+N R3R2R1R4+NR1 R2N+R1NR2+NR1X –−B FFFF−PFFFFF FCFFFSOO−NSOOCFFFFigure 1.1: Common Ions found in Ionic Liquids. Typical cation structures areshown across the top row and include the tetraalkylammonium cation, dialkylpyrro-lidinium cation, alkylpyridinium cation, and alkylimidazolium cation, where R1-4 in-dicate alkyl substituent groups. Typical anions are shown across the bottom row andinclude halides, tetrafluoroborate, hexafluorophosphate, and bis(trifluoromethane-sulfonyl)imideAside from being salts with relatively low melting points, there is nothing thatrestricts ILs to particular structural features or functional groups. ILs have earnedthe moniker “designer solvents” because the physical properties can be tailored bychanging the combination of ions, or by changing the connectivities or identities2of substituent groups within those ions. Similar tactics are used to create new ILsby introducing (or exaggerating) molecular asymmetry, creating size mismatchesbetween ions, delocalizing charges, elaborating side chains, and/or making fluoroussubstitutions.6Like traditional inorganic salts, ionic liquids are dominated by electrostaticinteractions, but because their melting points are significantly lower, they requireless thermal energy (and the specialized equipment designed to withstand high tem-peratures) to reach the liquid phase. With an accessible liquid phase and tunableproperties, ionic liquids are being investigated for use in a wide variety of applica-tions.7–9 Before giving a few examples of the potential applications of ILs, we discusssome of the aspects of finding, and optimizing, a task-specific IL.For a specific task, it is possible to prioritize desirable physical propertiesand performance attributes of the material. Task-specific ILs can integrate func-tion, design, and safety elements at the level of ion structure, which encourages aholistic approach to material design. However, the identification of promising ILcandidates with particular physical properties is not trivial.10–12 Furthermore, de-termining which IL is best (as opposed to well-suited) for a particular applicationis even more difficult.13Early in the screening process, IL candidates are eliminated for undesirableproperties such as high toxicities and corrosivities.14 Thermal stability (decomposi-tion) and reactivity (potential degradation) are considered for the operating condi-tions and expected life cycle of the application. The thermal conductivities and heatcapacities are also important.14 Tunable properties such as the densities, viscosities,hydrophobicities,15,16 electrical conductivities, liquid phase range, hygroscopic ten-dencies, surface and bulk structures, and/or solubilities are tailored into the desired3range for the particular applications. Wide electrochemical windows—the voltagerange limited by oxidation and reduction of the ions—are also a common featureand desirable attribute of ILs for electrochemical applications.The range of IL properties far exceeds those of conventional organic sol-vents.17 ILs have been called “green” solvents because of their low volatilities andthe environmentally attractive possibility of being able to recycle/reuse ILs withrelative ease, in comparison to more traditional volatile organic solvents. The desig-nation of ILs as environmentally friendly solvents does not address issues of potentialtoxicity or biodegradability.3 For solvent use, the bounds on the liquid range, vis-cosity, solubilities, and non-reactivity are the key tunable parameters. ILs providea reaction environment with more electrostatic character, which influences the solu-bilities of organic and inorganic materials. IL solvents open up synthetic and kineticpathways that may be more efficient (and diverse) than traditional organic solventchemistry. Reactivity in terms of reaction rates, yield, product selectivity, and evenenantioselectivity are all subject to change with different solvents.3,18,19 Chiral ILshave been used successfully in asymmetric synthesis.20ILs are also promising as electrolytes.14 The desirable properties of elec-trolytes include a wide electrochemical window, high conductivity, wide operatingtemperature range, and for safety reasons, low flammability and low volatility.14,21As the scale of battery-powered devices increases into vehicle-sized objects, the de-mand for safer batteries also increases.21 Large batteries based on strong electrolytes,such as H2SO4 (aq) and KOH(aq), would present safety hazards due to the amounts(and concentrations) involved. The safety measures required for containment aredue to extreme pH and reactivity, not necessarily the high conductivity. ILs arean enticing alternative. The tunable nature of IL properties can be exploited to4dampen reactivity but maintain sufficient conductivity for use as an electrolyte.ILs have been investigated for a variety of other novel applications. Wenote a small sample of them here. High thermal stabilities, densities, and heat ca-pacities also make ILs promising candidates as heat (energy) transfer and storagefluids.14,22–26 Long alkyl chains on ILs increase the likelihood of tail aggregation,making ILs attractive as surfactants.27,28 Tunable viscosities and high thermal sta-bilities make ILs attractive options for lubricants.29–31 ILs with low vapour pressurethat are capable of absorbing large amounts of CO2 are being investigated for car-bon sequestration.32–36 ILs can be tailored to have optimal solubilities for a targetedchemical species and (im)miscibility with traditional solvents, making ILs extremelyattractive as separation/extraction media.6,37–39 ILs can also serve as multifunctiontools in catalysis,33,40–44 where they can act as a solvent to bring together bothorganic and inorganic species, act as a catalyst or co-catalyst, or act as a ligandsource.42,44The peculiar physical properties of ILs are not yet well understood. Thecomplex relationship between molecular structure and physical properties is an im-portant research question from both theoretical and practical standpoints.10,45–47Structure-property relationships, based on the earliest known ILs, developed in arather empirical fashion.48 A goal of theoretical research is not only to explainthe trends in the physical properties that are observed experimentally, but also tounderstand the fundamental physics that governs IL behaviour and move towardpredicting the behaviour of potential ILs based on their chemical structures.15,49,50The Holy Grail for ILs would be being able to accurately predict all of the physicalproperties from the chemical structures of the ions alone.51.2 Computational Studies of ILsThe advances in computational resources are well timed for IL research.47Much simulation work is being done to find correlations between the molecular struc-ture of ILs and their physical properties13,51–55 and to determine accurate meltingpoints of model ILs.2,7,10,56Malvaldi and Chiappe studied structural and dynamical properties of a set ofthree salts with two-site cations and single-site anions.45 They found orientationalordering and structural heterogeneity in the liquid phase. Ganzenmu¨ller and Campconsidered a set of five salts with one- and two-site cations and single-site anionswith a size ratio of 1.3014:1.57 They located the critical points and resolved cationtranslational and rotational motion in the conductivity spectra. Roy and coworkershave developed and refined a more realistic, three-site coarse grain model for an alkylimidazolium cation and paired it with a single-site hexafluorophosphate anion.58,59Their studies of the liquid phase dynamics showed that the rotational diffusioncoefficient appears to decouple from the viscosity as the temperature decreases.Earlier simulations60,61 employing simple model ions have shown that sometrends in the dynamical properties of ILs can be correlated with a measure calledthe charge arm. The charge arm is defined by Kobrak and Sandalow as a mea-sure of the separation of the center of charge and the center of mass.62 It is oneof the quantities we will discuss in more detail in Chapter 2.4. The simulationstudies60,61 showed that ions with small charge arms display decreased viscosity andincreased diffusion and electrical conductivity, whereas ions with large charge armsshow a significant increase in viscosity, constant diffusion rates, and the electricalconductivity approaches zero. While we employ simple model ions with fixed charge6arms, the charge arm is not a fixed quantity in molecular ions with conformationalflexibility. Urahata and Ribeiro studied the heterogeneous dynamics of 1-butyl-3-methylimidazolium chloride using an atomistic model in MD simulations.63 Theyfound that the instantaneous charge arm acts as a structural signature of hetero-geneous dynamics with larger charge arms being associated with slower dynamicalbehaviour.631.3 The Melting Points of Ionic LiquidsThe melting point is a key physical property for IL applications and is thedefining attribute that separates ILs from molten salts.1 The fluid-solid transitionhas not received as much research attention as fluid-fluid transitions.64 The normalmelting point is the physical property we focus on in this thesis. Understandinghow a variety of structural modifications influence the melting point is a key stepin developing predictive tools for the rational design of ILs.12,47,53,55,65It is not well understood why ILs have such low melting points compared tosimple inorganic salts.66,67 The low melting temperatures are usually explained bythe differences in electrostatic interactions, namely the ion sizes, charges, and chargedistribution.66,68 There are many additional factors that nuance the relationship be-tween the Coulombic interactions and melting points such as competing interactiontypes and the associated length scales, geometric considerations, and conformationalflexibility.51,55,68 While it is difficult to unravel and isolate the contributions of a sin-gle factor,69 we outline a few of the generic features and the arguments for how theycontribute to the low melting temperatures of ionic liquids.7At the level of ion structure, molecular asymmetry is usually one of thefirst reasons offered to explain the low melting points.65 Asymmetric cations aremore prevalent than asymmetric anions, although the anions are becoming moreexotic. The ion asymmetry introduces an orientational dependence that can leadto reorientational degrees of freedom. Ion asymmetry may prevent the electrostaticordering (lattice packing) that drives crystallization, thereby lowering the meltingtemperature.7,54,63,70 The work presented in Chapters 3, 4, and 5 of this thesis studyhow redistributing the cation charge affects the melting temperature.Ionic size is another clear factor that affects the melting temperature, withlarger ions having weaker electrostatic interactions (electrostatic energies ∝ 1/r)and therefore lower melting temperatures.64,71 The work presented in Chapter 4 ofthis thesis examines the impacts of ion size and ion size ratios, and also examinesthe compounding effects of ion size, ion size ratios, and cation charge redistribution.Conformational flexibility is another reason put forward to explain the lowmelting points. The mix of conformers with different shapes and charge distributionsmay frustrate the crystal packing.7,52Combining a cation with an anion generates more factors to consider. Amismatch in ion size affects the packing and freezing behaviour. Large ions aremore likely to have significant van der Waals interactions than smaller ions. The sizedifference can lead to a prefreezing or premelting phenomenon, where the diffusivemotion of the smaller anions persists at lower temperatures than that of the largercations. We explore the impact of size ratio in Chapter 4, and find salts with fast ionconductor or plastic crystal phases. Experimentally, salts with plastic crystal72–74and fast ion conductor75 phases of ILs have been studied for potential use as solidstate electrolytes.8Monoatomic ions are spherical, whereas molecular ions are typically aspheri-cal. Just as combining two ions of different sizes introduces size ratio, combining twoions with distinct shapes introduces the possibility for shape “mismatches.” Shapemismatch between the ions is often cited as a reason for the low melting points ofILs. If the mutual fit of the ion shapes is poor, the mismatch can limit the packingoptions.20 Even if the ions find a densely packed structure that accommodates bothion shapes, the local arrangement of ions may lead to unfavourable interactions thatdestabilize the crystal structure.Another feature that is introduced once we consider the cation/anion combi-nation is the possible formation of ion pairs. Ion pairing has been correlated to lowglass transition temperatures.48 We also see strong correlations between directionalion pairs and low glass transition/melting temperatures in Chapters 3, 4, and 5. Ionsthat form strong, long-lived directional ion pairs may not interact very strongly withother neighbouring ions. The thermodynamic driving force behind solidification isthe lower enthalpy in the solid phase. The lower enthalpy usually depends on par-ticles creating strong, favourable interactions with multiple near neighbours, whichis unlikely to be the case for salts that form strong, directional ion pairs.Turning our focus to the solid, the bulk properties of the solid phase areanother complicating factor for ILs. The crystal structure itself must be known inorder to find a relationship between melting temperature and ion structure.76 Thestrength (large lattice energy), or rather lack thereof, of the crystal structure is of-ten cited as a reason for low melting ILs.51 ILs, with their diverse properties, canvitrify,15,48,77 crystallize,51 crystallize with orientational disorder, crystallize as plas-tic crystals,73,78–80 crystallize as fast ion conductors,75,81 form liquid crystals,77,82,83and/or exhibit polymorphism.15,55,77,84 We have observed all of these solidification9behaviours with our simple salt models except for the formation of liquid crystals.We would have to move away from spherical ions to induce liquid crystal behaviour.The dynamics in the solid phases of ionic liquids, including oscillations, reorien-tations, rotations, conformational changes, are especially active near the meltingtransition. The rich diversity of thermal motions in the solid add nuance to locatingthe onset and end of the solid-liquid (s-l) transition.It is important to acknowledge the interdependencies among the levels of ionstructure, pair energies, bulk liquid structure, and crystal packing that make it verydifficult to isolate qualitative trends in melting temperatures.52 The factors givenhere, and their interdependencies, support the idea that the morphologies of ILsare more complex than that of molten salts,1 and particularly, the origin of the lowmelting temperatures is an important research question.1.4 Experimental Melting PointsNormal melting points of pure ILs are a challenge to determine accurately.The main difficulty lies in obtaining pure, bulk samples (free of other ions, water,and carbon dioxide gas) and using a method that is reliable and accurate. IL samplesare commonly dried under vacuum for hours to remove water and then stored inan inert atmosphere before measuring the melting point. It has been reported thatdissolved carbon dioxide gas in tetrabutylammonium tetrafluoroborate can reducethe melting point by 120 K, from 429 K to 309 K at 150 bar of CO2.26,85Experimental melting points are subject to the compounding issues of samplepurity, variations in sample preparation procedures, and the accuracy of the methodemployed. Differential scanning calorimetry and visual observation are two common10approaches for determining melting points. Reported experimental melting pointsfor the same salt measured using different methods (and samples) do not necessarilyagree within the given uncertainties of the method and are commonly separated by25-35 K.26 Reliable melting point data is limited to the popular ILs that havebeen characterized and studied extensively. One common ionic liquid, 1-hexyl-3-methylimidazolium bis[trifluoromethyl)sulfonyl]amide, transitions to the liquid statefrom the liquid-quenched glass around 265 K, and from the crystal around 271K.26,86 IUPAC has proposed to use this common ionic liquid as a reference IL withits relatively narrow melting point range of 266-272 K.26,87 The limited availabilityof accurate melting point data makes establishing empirical trends difficult.The discrete nature by which molecular ions can be modified means that,within a homologous series of salts, substituting one (set of) neutral atom(s) foranother can change the ion geometry (size and shape), mass and mass distribution,and the charge distribution of a molecular ion. Isolating how a single physicalattribute affects the melting point using a series of salts with molecular ions andexperimental data increases the complexity of the task dramatically.In order to connect the results obtained here for coarse grain model salts toexperimental ones, we take a brief look at the melting point trends across a seriesof alkyl-substituted ammonium bromide salts in Chapter 6.111.5 Simulation Methods for Estimating MeltingTemperaturesEquilibrium melting temperatures are well defined thermodynamically. Thechemical potentials of both phases vary slowly with pressure. Determining theprecise coexistence points using computer simulations is challenging, especially formolecular crystals which are complex.88The s-l transition is impossible to characterize completely without knowl-edge of the solid phase structure. In this work, most of the salts (about 84%)spontaneously crystallize when cooled from the liquid phase, while a minority ofsalts become trapped in “glassy” states. We explore the solid state structures ofthe glassy salts by preparing and heating common crystal structures using NPTmolecular dynamics. Another method of determining crystal structures was alsodeveloped for this project. A genetic algorithm was written to find low energy unitcells of salts. The calculation details are described in Appendix A and the programis described in Appendix C.Different approaches have been developed to estimate melting points fromsimulations and, not surprisingly, the more accurate methods tend to be computa-tionally intensive.69,89,90 A collection of methods has been recently reviewed in anexcellent work by Zhang and Maginn.89 The two methods we use, the hysteresismethod and two-phase simulations, are described briefly here and in more detail inChapter 2.We use the hysteresis method,91 as an initial estimate of the melting point inthis work. The hysteresis method determines the melting point from the transition12temperatures of freezing the liquid and melting the solid. The freezing and meltingtransitions do not occur at the same temperature due to hysteresis. Hysteresisdescribes the phenomenon of a lag in response to a change in conditions or forcesacting upon a system. The state of a hysteretic system depends on the historyof the system as well as the current conditions. When the liquid is cooled, thesystem will remain liquid below the equilibrium freezing point (supercooling) andwhen the solid is heated it will remain solid above the equilibrium melting point(superheating). Hysteresis is a non-equilibrium phenomenon, and the dominantcause of the hysteretic behaviour around the s-l transition is the free energy barrierassociated with nucleation.91One disadvantage of using the hysteresis method for the ionic systems sim-ulated in this thesis is that the thermal width of the hysteresis is large and spanshundreds of degrees. To refine the melting temperature estimate, we use NPT MDsimulations that start as a composite of two phases. The final solid and liquid con-figurations from the hysteresis simulations are used to create composite supercellswith two s-l interfaces. We use the two-phase simulations to narrow the transitionrange to 50 K.In this work, we are interested in the melting point trends as the ion param-eters are systematically varied and not the exact melting points of individual salts.Therefore, our estimates of the melting points of individual salts are not the mostaccurate possible. However, the methods employed are expedient and sufficient foridentifying trends across sets of salts.13Chapter 2Models and Methods ∗We employ molecular dynamics (MD) simulations to estimate the meltingpoints of sets of coarse grain model salts. In this chapter, we first present thedetails of the ion models, and then the simulation methods.2.1 Ion TopologiesWe use a coarse grain approach in this work because it provides insightsinto the influence of particular molecular features. By reducing the complexityof the ions, relationships among model parameters and observables, particularlymelting points, become more prominent. Coarse grain models have been used tostudy the fluid phases of ILs and have been successful at reproducing some observedbehaviour.60,61,92∗Parts of this chapter have been published. E. K. Lindenberg and G. N. Patey, “How distributedcharge reduces the melting points of model ionic salts,” J. Chem. Phys. 140, 104504 (2014).14All of the salts we consider can be viewed as extensions of the so-calledrestricted primitive model (RPM). The RPM is an equimolar mixture of equal-sized, charged hard spheres, and is perhaps the simplest model of inorganic salts.93We modify this model by introducing Lennard-Jones (LJ) interactions, by changingthe cation-anion size ratio, and by redistributing a fraction of the cation chargeaway from the ion center. Charge redistribution is not an applicable extension ofthe RPM for monoatomic ions, where the center of mass and the center of chargenecessarily coincide, however, it is applicable for molecular ions, which are commonin real ILs.Each ion is monovalent and has a single LJ interaction site that, for ourpurposes, defines the ion center. The LJ parameters are specified in Chapter 2.2.Each ion has a total mass of 1.993×10−25 kg (120 amu), and the mass distributionsare discussed further in Chapter 2.3.We present a brief overview of each ion geometry before specifying the com-plete ion parameters. A summary table is included at the end of this section wherewe list the salts studied in the following chapters.The 1C IonThe 1C ion is closely related to the ions in the RPM. The ions have a singleinteraction site located at the ion center that carries the entire unit charge. Unlikethe RPM ions, the 1C ions also interact via a LJ potential. The 1C - 1C salt actsas a reference salt for our work and is included in each chapter. The deviations inthe melting points of the other salts are measured relative to the 1C - 1C salt.15The 2L IonIons with some charge located at an off-center site are labelled 2L, indicatingtwo interaction sites in a linear arrangement, as shown in Fig. 2.1.dFigure 2.1: A diagram of the 2L model ions. Each ion carries a total charge of ±1e,where e is an elementary charge. The site at the ion center interacts via a Lennard-Jones potential and, in most cases, an electrostatic potential. The off-center siteinteracts via an electrostatic potential only. The two parameters for the 2L ionsare the amount of charge at each interaction site and the distance between the twosites, d.The 2L ion parameters are the distance between the sites d and the fractionof charge on each site. We examine 2L ions with four different charge distributions,with 1/3, 1/2, 2/3, and all of the unit charge on the off-center site. The ion labelsare 2Lfc-d, where fc is the fractional charge on the off-center site as a roundedpercentage, and d is the separation distance between the two interaction sites onthe same ion, expressed in units of 0.01 nm. For example, a 2L67-12 ion has 2/3of the total ion charge on the off-center site, which is located 0.12 nm from the ioncenter.16The 3s IonsThe 3s ions are a set of ions with three interaction sites. The third ion siteis another off-center charge site, which, like the off-center charge site in the 2L ions,is located a distance d from the ion center. The third site introduces another modelparameter, θ, the angle made between vectors from the ion center to the two off-center sites. We systematically vary θ across seven values: 60◦, 90◦, 105◦, 120◦, 135◦,150◦, and 180◦. An example of a labelled cation is shown in Fig. 2.2, and minimalistdiagrams of the entire 3s ion set are shown in Fig. 2.3.CMCC1CC2dd θ = 120◦Figure 2.2: A fully labelled schematic of the 3s(120◦) model ion.The 3s ion geometry labelling follows a convention similar to that of the 2Lions, given as 3s(θ)fc-d, which specifies the ion geometry, fractional charge distri-bution, followed by the extent of the charge separation d. The angle θ is in degrees,fc denotes the percentage of charge on interaction site CC1, and d follows the sameconvention as the 2L ions. An example of a full 3s ion label is 3s(105◦)33-18. Theion label gives the number of sites (3s), the angle between the two off-center sites(105◦), the amount of charge at CC1 [fc = 33, q(CC1) = 1/3], and finally, thedisplacement distance d is given in units of 0.01 nm (0.18 nm).17In the 3s - 1C salts we consider, the off-center interaction sites are both thesame distance d from the ion center. The charge distribution fc specifies the amountof charge on one off-center site (CC1), while the other off-center site (CC2) alwayscarries 1/3 of the unit charge. In Chapter 5, we study two variations in chargedistribution. In the first case, the CC1 site carries 1/3 of the charge (fc = 33) andthe ion center carries 1/3 of the unit charge. In the second case, CC1 carries 2/3of the unit charge (fc = 67) and the ion center is uncharged (but still interacts viathe LJ potential).183s(60◦) 3s(90◦)3s(105◦) 3s(120◦) 3s(135◦)3s(150◦) 3s(180◦)Figure 2.3: Diagrams and labels of the 3s ions. The cationic version of these sevenions, each paired with the 1C anion, are studied in Chapter 5.19Ion Size RatiosWe have defined the general features of the ion topologies in the precedingsections, and now turn to the combination of two ions to create a salt. The rela-tive size ratio of the two ions, as well as the absolute ion sizes are both relevantparameters to be varied. With the exception of the salts considered in Chapter 4,the cations and anions are the same size, with a LJ diameter σ = 0.50 nm.CM CC AMdσ±σ+σ−Figure 2.4: A schematic of the size-asymmetric model salts. The interaction sites,CM, CC, and AM denote the cation center, cation off-center charge site, and anioncenter, respectively. The size ratios are varied in Set A by fixing σ± to 0.50 nmwhile varying σ+ and σ−. In Set B, the size ratios are varied without a constrainton σ±; instead σ− is held constant while σ+ and σ± increase. In addition to varyingthe size ratio, we also examine two cation charge distributions where either all ofthe unit charge is located at CM, or 2/3 of the unit charge is located off-center atthe CC site, at a distance d from the ion center.20The cation:anion radius ratio is controlled by the length parameter of the LJinteraction, σ. Allowing the ion sizes and resulting size ratios to vary creates threedifferent length scales to consider, σ+, σ−, and σ±, which affect the cation-cation,anion-anion, and cation-anion interactions, respectively. We vary the cation:anionradius ratio in two different ways, generating two sets of size-asymmetric salts. InSet A, the distance between cation-anion centers σ± is held fixed, which isolates theinfluence of size asymmetry and minimizes the changes in strength of the attractiveelectrostatic interactions. In Set B, σ+ is increased (which also increases σ±) whileσ− is held constant. The salts in Set B are somewhat more practical representationsof ILs because σ− = 0.50 nm was initially chosen to be consistent with the “size”of typical IL anions. We explore seven different size ratios, ranging from 1:1 to 3:1.The parameters for each size set are given in Table 2.1. Note that the size-symmetriccases, A100 and B100, are identical; we refer to these salts using the A100 label only.Table 2.1: The two sets of size parameters of the salts considered in this work.In Set A, σ± is constrained to 0.50 nm, while in Set B, the anion radius is heldconstant at 0.50 nm and σ± is unconstrained.Size Label Set A Set BRatio σ+ (nm) σ− (nm) σ+ (nm) σ− (nm)1:1 100 0.500 0.500 0.500 0.5001.33:1 133 0.571 0.429 0.667 0.5001.67:1 167 0.625 0.375 0.833 0.5002:1 200 0.667 0.333 1.000 0.5002.33:1 233 0.700 0.300 1.167 0.5002.67:1 267 0.727 0.272 1.333 0.5003:1 300 0.750 0.250 1.500 0.50021In Chapter 3, we consider five different sets of charge distributions in size-symmetric salts. In Chapter 4, we vary the size ratios and focus on two cation chargedistributions, where the entire unit charge is located at the cation center (labelled1C), or where 2/3 of the unit charge is displaced from the cation center (labelled2L67). We chose the 2L67 charge distribution because the size-symmetric 2L67 - 1Csalt series was one of two that obtained melting points of 50% of the size-symmetric1C - 1C salt in Chapter 3. Also, we found that the series of size-symmetric 2L67- 1C salts captured the typical behaviour of the salts we observed over the largerset of charge distributions. We vary the charge distribution on the 2L67 cations bychanging d, which ranges from 0.06 nm from the cation center up to the edge.The labels for the salts give the size ratio information first, followed by thecation and anion geometries. For example, the A233 2L67-18 - 1C model salt indi-cates the size ratio as A233, meaning the salt is in Set A, with a size ratio of 2.33:1,(the number gives the cation size relative to the anion as a percentage). The cationgeometry is specified next as 2L67-18, meaning 67% of the unit charge is displacedd = 0.18 nm from the cation center. Finally, the last term in the label specifies theanion geometry as 1C. The vast majority of the salts considered here are combina-tions of a 2L67 cation with a 1C anion; it is the size ratio information and d thatare the key variables in Chapter 4. Note that the size ratio (A100) is excluded fromthe size symmetric salt labels in Chapters 3 and 5.22As the ion sizes change, the cutoff distance for the short-range interactionsalso changes. The short-range interaction cutoff distance is maintained at approx-imately 3.2σ+ and increases as the cation size increases. The cutoff distances aregiven in Table 2.2.Table 2.2: The short-range interaction cutoff values (Rcut) for the size-asymmetricsalts. The LJ diameter of the cation, σ+ is larger than the anion for all of thesize-asymmetric salts and was used to determine Rcut as Rcut ≈ 3.2σ+.Size Ratio σ+ RcutLabel (nm) (nm)A100 0.500 1.70A133 0.571 1.83A167 0.625 2.00A200 0.667 2.14A233 0.700 2.24A267 0.727 2.33A300 0.750 2.40B133 0.667 1.83B167 0.833 2.00B200 1.000 2.14B233 1.167 2.24B267 1.333 2.33B300 1.500 2.40232.2 Ion Interaction ParametersThe ions have a combination of electrostatic and LJ interactions because wewant to include both ionic and molecular character in our model ions. The ionsinteract through the pair potentialu(rij) = 4εij[(σijrij)12−(σijrij)6]+ 14pi0∑α , βqiα qjβriαjβ, (2.1)where εij and σij are the LJ energy and length parameters for the LJ interactionbetween ions i and j, rij is the distance between the LJ sites on ions i and j, 0 isthe permittivity of free space, qiα is the charge on ion i at site α. The riαjβ valueis the distance between site α on ion i and site β on ion j, where both sites carrycharge.Each ion has one LJ interaction site at its center. The LJ energy parameteris identical across all models; εij = 6.0× 10−21 J. The value of εij is consistent withearlier work60,61,92 and is roughly consistent with the model parameters of a typicalhydrocarbon, neopentane.94The ion topologies were set up with each ion site as an “atom” with a label,mass, charge, and LJ parameters. The connectivities between interion sites were setup as distance constraints. An example topology file for the A100 2L100-16 - 1Csalt is given in Appendix F.In Chapters 3 and 5, the LJ length parameter is held constant across allion models, σ = 0.50 nm. The ion diameter of 0.50 nm approximates the ion“size” in some common ILs, such as the PF –6 or BF –4 anions. In Chapter 4, theLJ length parameter changes to create salts with cation-anion size ratios varying24from 1:1 to 3:1. The cation-anion length parameter, σij, is calculated using theLorentz-Berthelot mixing rules, σij = (σii + σjj)/2.In the charge-centered case, the pair potential is spherically symmetric anddepends only on the distance from the ion center. When d increases in the 2L ions,part of the charge is moved closer to the ion surface. For the 2L - 1C salts, thelocation and depth of the potential minima changes dramatically from the 1C - 1Ccase. Three dimensional (3D) surface plots of the pair potential for select 2L67 - 1Csalts are shown in Fig. 2.5. Each plot shows the potential energy surface around aparticular cation, as experienced by a 1C anion.251C - 1C 2L67-06 - 1C−500.0−400.0−300.0−200.0−100.00.0pairpotential(kJ/mol)2L67-10 - 1C 2L67-14 - 1C2L67-18 - 1C 2L67-22 - 1CFigure 2.5: 3D potential energy surfaces around the cation for the 1C anion in aseries of size-symmetric ion pairs. Each cation center is shown at the origin withthe circle depicting the repulsive core (where u(rij) > 0). The off-center charge siteon the 2L67 cations is aligned with the +x axis.262.3 Ion MassesFor the 2L and 3s cations, the ion center defined by the LJ site is distinctfrom the ion’s center of mass because a fraction of the mass is located at the off-center site(s). The amount of mass at each interaction site depends on the cationdisplacement distance, d. The distribution was selected to minimize fluctuations inthe moments of inertia as d varies. The distribution of mass affects the dynamics,but does not affect equilibrium properties such as the normal melting points. Themasses assigned to the ion center (labelled CM) and each off-center ion site (labelledCC) are given in Table 2.3.Table 2.3: The mass distribution for the ions with off-center interaction sites.The total mass of each ion is 120 amu (1.993×10−25 kg), and the mass distributionover the interaction sites varies with the separation distance d. For the 3s ions,both off-center sites were assigned the same mass.2L Ions 3s Ionsd m(CM) m(CC) m(CM) m(CC)(0.01 nm) (amu) (amu) (amu) (amu)0 120.00 — 120.00 —6 70.00 50.0010 102.00 18.0014 112.00 8.0018 115.00 5.00 115.00 2.5022 116.75 3.2526 117.75 2.2530 118.25 1.7534 118.70 1.30≥ 42 118.70 1.30272.4 Charge ArmThe charge arm, defined by Kobrak and Sandalow62 as |qoff center| d, is a usefulcombination of both parameters that occur for 2L ions. The notion of charge armis a defining feature in the set of 2L - 1C ion models considered here. In the presentwork, we use the normalized charge armLc =|qoff center| d|qion|Rion, (2.2)where qion is the total ion charge and Rion is the radius of the ion. In this work, allthe ions considered are monovalent, so |qion| = 1e. With the exception of the workin Chapter 4, Rion = 0.25 nm. In Chapter 4, Rion ranges from 0.25 nm to 0.75 nm.The charge arm characterizes the geometry of one ion only, and does not considerchanges to the other ion. In Set A of Chapter 4, the cation charge geometry andboth the cation and anion radii change, so a measure that reflects only the change inthe cation may be misleading. Rather than using Lc for Set A, we simply correlatethe melting temperatures with d for each size ratio (Fig. 4.5). In Set B of Chapter 4,however, the anion size is fixed, and it is instructive to consider the changes to thesalt through the normalized charge arm (Fig. 4.7).A model cation with a smaller (larger) charge arm is somewhat analogousto a heterocyclic cation with a shorter (longer) alkyl side chain. In the model ions,displacing the charge away from the ion center creates a stronger electrostatic inter-action on one side of the ion and effectively “exposes” the van der Waals interactionon the opposite side.282.5 Minimum Potential Energies of Isolated IonPairs, U IPisoThe minimum potential energy for an isolated pair of ions, U IPiso , is anotheruseful quantity for characterizing the salt models. The pair potential considers thegeometry of both ions in the salt, rather than just one, which permits comparisonacross salts where the geometries of both ions change. A linear arrangement ofCM-CC-AM, as depicted in Fig. 2.4, is used to compute U IPiso .Once we have determined a crystal structure for a salt, we use the ratio of theaverage configurational energy per ion pair in the crystal at the melting temperatureto the minimum potential energy of an isolated ion pair, U¯ IPs /U IPiso, to provide someinsight into what fraction of the crystal energy is coming from ion pairing rather thanfrom the packing arrangement. This ratio is somewhat analogous to the Madelungconstant for point charges on a perfect lattice. The Madelung constant is a factorthat relates the total electrostatic potential of a specific ionic crystal structure tothe minimum distance between two counterions, r0, byUES = Mq+q−4pi0r0, (2.3)where M is the Madelung constant. By rearranging Eq. (2.3) and substitutingin the minimum potential energy of an isolated ion pair, U IPiso = q+q−/4pi0r0, theexpression for the Madelung constant becomesUES/U IPiso = M. (2.4)29Like the Madelung constant, U¯ IPs /U IPiso is a measure of how much stability isgained in the solid phase relative to a collection of isolated ion pairs.65 At a ratio of1.0, the solid would have no energetic advantage over an equal number of isolatedion pairs. However, U¯ IPs /U IPiso is unlike the Madelung constant in that it includes theLJ interactions, thermal fluctuations, and any lattice defects present.2.6 Outline of the Salts Studied in This ThesisIn Chapter 3, we examine a set of 2L - 1C salts made by combining a 2Lcation with a 1C anion. The size-symmetric 1C - 1C salt, like the RPM, is perfectlysymmetric with respect to sign reversal on the charges. The 2L - 1C salts extendthe simple models by introducing asymmetric charge distributions, which are morerealistic for modelling ionic liquids. The choice of redistributing the cation charge isarbitrary and reversing the charge signs would produce a 1C - 2L salt with identicalphysical properties. Chapter 3 examines the effect of varying the charge distributionon the 2L cations. The ions are the same size and the anion is always a 1C ion.In Chapter 4, we examine how the absolute ion sizes and size ratios affectthe melting temperatures of 2L67 - 1C and 1C - 1C salts. The coarse grain modelsalts considered in Chapter 4 mimic two features commonly found in ionic liquids,size asymmetry and distributed cation charge.In Chapter 5, we study the set of 3s - 1C salts to further explore the affectof cation charge distribution on the s-l transition. The ions are size-symmetric andalways paired with a 1C anion. While we vary θ in this chapter, the displacementdistance d is held fixed at 18 (0.18 nm). The displacement distance is far enough tohave an effect on the melting temperature, yet not so far that ion pairing dominates30the ion interactions.Table 2.4: Table of the salts studied in each chapter of this thesis. Each chapteralso includes 1C - 1C salt(s) as the reference salt(s).Chapter Cation(s) Anion(s) Ion Sizes d valuesσ (nm) (0.01 nm)3 2L 1C 0.50 0 - 244 (Set A) 2L67 1C 0.25-0.75 0 - 264 (Set B) 2L67 1C 0.50-1.50 0 - 585 3s 1C 0.50 182.7 Hysteresis Method for Estimating MeltingTemperaturesSince we are interested in melting point trends as the ion parameters aresystematically varied rather than determining the exact melting point of individ-ual salts, our estimates of the melting points are not the most accurate possible.However, the methods we employ are expedient and sufficient for identifying trendsacross sets of salts.One relatively simple method for estimating melting points is the so-calledhysteresis method proposed by Luo et al.91 A hysteresis loop is evaluated by su-percooling the liquid until it freezes at a temperature below the thermodynamicmelting point and superheating the solid until it melts at a temperature above thethermodynamic melting point. The melting point, Tm, is estimated using the limits31of supercooling and superheating as91Tm = T+ + T− −√T+T− , (2.5)where T+ is the highest temperature of the superheated solid and T− is the lowesttemperature of the supercooled liquid.When a system is supercooled or superheated, the difference in chemicalpotential provides thermodynamic incentive to freeze or melt. Under isobaric con-ditions, the rate of heating or cooling determines the maximum extent of the super-heating or supercooling, and has been established as91,95β ∝ θc(1− θc)2 (2.6)where β is a dimensionless nucleation barrier parameter and θc is either the maxi-mum extent of superheating (θ+ = T+/Tm − 1) or supercooling (θ− = 1− T−/Tm).The nucleation barrier parameter, β, should be identical for both melting and freez-ing, soθ+(1− θ+)2 = θ−(1− θ−)2. (2.7)From Eq. (2.7), an expression for the melting point, Tm, [Eq. (2.5)] is obtained byusing the quadratic formula.The hysteresis method has been used to obtain reasonably accurate estimatesof the melting points of a pure LJ system91 and of molecular solids.88,96 Since weare using rather large temperature increments of 50 K, we employ Eq. (2.5) as aninitial estimate.32We started each simulation by equilibrating the salt in the liquid phase. Theinitial simulation temperature was 1600 K for most salts. Some salts evaporated at1600 K, so the initial temperature was lowered until we found a stable liquid. Themolten salts were then cooled using serial NPT runs with temperature differencesof 50 K. The freezing of the supercooled liquid was identified by significant decreasesin the average configurational energy and volume, and confirmed by inspection ofthe radial distribution functions and mean square displacements. The solid wasthen heated following an analogous procedure. For most of the salts, the solid-liquid transition temperatures from the cooling and heating cycles show evidence ofhysteresis.A few of the salts (about 16%) did not crystallize under our simulation condi-tions and time scales, but tended to form “glassy” states as they were cooled. Thesesalts are discussed in the relevant results chapters. The vitrification creates threerelated issues if one wishes to estimate Tm using Eq. (2.5). We cannot determine thevalues of T− or T+ directly from cooling and heating curves, and we do not knowwhich crystal structure(s) exist below the glassy state. Zheng et al.88 have arguedthat a reasonably good estimate of T− is given by a glass transition temperature,Tg, defined as the temperature where the diffusion coefficient falls below 10−11 m2/son the cooling curve, and we follow this prescription here. The value of T+ can bedetermined by heating a prepared crystal until it melts, if the appropriate crystalstructure is known.To find candidate crystal structures, we set up common crystal types (suchas CsCl, NaCl, NiAs, and wurtzite) for each glassy salt. The initial crystal densitywas set at ρ∗ = 1.1, where ρ∗ = Nσ3/V (N is the number of ions, σ is the LJdiameter, and V is the volume of the simulation cell). The configurational energies33were minimized, and the crystals were heated from T = 50 K to determine T+. Uponheating, some of the solids rearranged into more stable crystal structures, which arediscussed in the relevant results chapters. The configurational energies of the newcrystal structures were minimized for each of the glassy salts and the crystals wereheated beginning at 50 K, as with the other prepared crystals.2.8 Two-Phase Simulation Method forEstimating Melting TemperaturesTo refine melting point estimates, we simulated supercells of N ≈ 1600(Chapter 3) or N ≈ 4000 (Chapters 4 and 5) particles prepared with two solid-liquid interfaces. The two-phase simulation setup, shown in Fig. 2.6, is expectedto eliminate some of the kinetic barriers to homogeneous nucleation, and narrowthe range of our melting temperature estimates with respect to the thermal rangeof the hysteresis loops (T− to T+). From the hysteresis simulations, we have liquidand solid configurations in 50 K increments between T− and T+. We construct as-l-l-s composite from the solid and liquid configurations for each simulation tem-perature between T− and T+. We use each composite as the initial configuration fora short NV T equilibration simulation and then a longer NPT simulation. Duringthe simulations, each system usually evolves into a single phase.The series of independent two-phase simulations carried out in 50 K incre-ments between T− and T+ determine another two transition-related temperatures.The highest temperature at which a two-phase simulation evolved into a solid isdenoted Ts, and the lowest temperature at which a two-phase simulation evolvedinto a liquid is denoted Tl. The difference between Ts and Tl is usually 50 K, the34temperature grid spacing between independent two-phase simulation runs. In a fewcases, the two-phase simulations did not evolve into a single-phase; both phases per-sisted through the entire simulation. The persistence of both phases indicates thatthe salt is at or near its melting temperature, rather than clearly above or belowit. For these exceptions where a two-phase simulation ended apparent coexistence,we take the apparent coexistence temperature as Ts and the first temperature thatgives a pure liquid as Tl.An example configuration of the two-phase simulations is shown in Fig. 2.6.The central simulation cell and the simulation cell with periodic images in twodirections are both shown. The boundaries around the central simulation cell aredrawn in black. The periodic images are shown to emphasize that the number ofparticles and their arrangement in the s-l-l-s composite cell should emulate bulkbehaviour for each phase.35Figure 2.6: Snapshots of the initial configuration for the two-phase simulations forthe size-symmetric 1C - 1C salt at T = 1250 K. The full simulation cell (top) and thesimulation cell with periodic images in two directions (bottom) are shown. Periodicboundary conditions were used for all of the molecular dynamics runs in this thesis,including the two-phase simulations. The cations and anions are shown in blue andred, respectively. The ions are not shown to scale to highlight the arrangement ofthe crystal.362.9 Simulation DetailsThe hysteresis loop simulations were done using Gromacs 4.5.497–99 in Chap-ters 3 and parts of 4, whereas Gromacs 4.6.2 was used for some of Chapter 4 and allof the simulations in Chapters 5. The change in Gromacs versions coincided with asystem upgrade on Compute Canada’s computing cluster Orcinus.To start the single phase simulations in the hysteresis loops, the coarse grainions were arranged in a CsCl crystal structure at a reduced density of 0.8 and equi-librated at sufficiently high temperatures to maintain the liquid phase. The systemswere cooled and heated by taking the final frame of one 4.0 ns NPT simulation at agiven temperature and using it as the initial configuration for the next 4.0 ns NPTsimulation at a temperature 50 K lower or higher. Example Gromacs input files(*.top, *.mdp, *.gro) and the execution script (*.pbs) are included in Appendix F.Each simulation in the hysteresis cycles was executed in the NPT ensemblewith a pressure of 1.0 bar. The number of particles in each simulation varies. InChapter 3 we used N = 432 particles (216 ion pairs) whereas in Chapters 4 and5, we used N = 1024 particles (512 ion pairs). We did not find any significantdiscrepancies in our final results based on the version of Gromacs used, or anysignificant system size dependencies.All simulations were done using custom forcefields to specify the ion topolo-gies as described earlier. A leap-frog integrator propagated the system through timewith a time step of 1.0 fs for all runs. The trajectories were 4.0 ns long with thefirst 2.0 ns for equilibration. Structural and dynamic properties were analyzed overthe final 2.0 ns, except for radial distribution functions (rdfs), which were analyzedover the last 1.0 ns.37We used a Nose´-Hoover thermostat100–102 and a Parrinello-Rahman baro-stat103,104 throughout. The Nose´-Hoover thermostat regulated the temperature ofthe system using two coupling groups, one for the cations and one for the anions.The relaxation time of the thermostat was 0.10 ps. An isotropic Parrinello-Rahmanbarostat was used for the runs in cooling and heating cycles with a relaxation timeof 5.0 ps; the compressibility of the salt was set to 3.0×10−5 bar−1 along the threeprimary axes, and 0.0 otherwise.Periodic boundary conditions and the minimum image convention were usefor the short-ranged interactions, which were spherically truncated at 1.70 nm for thesize symmetric salts, and at approximately 3.2σ+ for the size asymmetric salts. Thecutoff radii used in Chapter 4 are specified in Table 2.2. Long-range corrections wereapplied to the LJ energy and pressure. The long-ranged electrostatic interactionswere accounted for using the Particle Mesh Ewald (PME) method105 with a Fourierspacing of 0.20 nm. The constraints on “all-bonds” between interion sites weremaintained using the Lincs algorithm.106The simulation procedure for heating the prepared crystal structures is thesame as the cooling and heating procedure, but uses a larger number of ions (N =686 to 2240) and an anisotropic Parrinello-Rahman barostat to maintain the pressurein the solid phase.While the hysteresis simulations were split across two versions of Gromacs,the two-phase simulations were done with Gromacs 4.6.2, with the exception of thework in Chapter 3. There is some duplication of salts across chapters, particularlythe size-symmetric 2L67 - 1C salts, which were studied in both Chapters 3 and4. For the duplicated salts, the simulations were carried out within each set forconsistency. The results from Gromacs 4.6.2 with N = 1024 ions were consistent38with the earlier results from Gromacs 4.5.4 with N = 432 ions for all salts.We used the final configurations from the hysteresis simulations to create asupercell for the two-phase simulations at each temperature between T− and T+.The liquid configurations were cropped in two directions to match the dimensionsof the simulation cell of the solid, and excess ions were removed as necessary topreserve electroneutrality. A small space (equivalent to σ+) was left at each s-linterface. The configurational energies were minimized to remove poor interactionsacross the new interfaces, and the system was re-equilibrated in the NV T ensemblefor 0.25 ns before running under NPT conditions for 4.0 ns. A run length of 4.0 nswas sufficient for most systems to evolve into a single solid or liquid phase. Periodicboundary conditions were also used for the two-phase simulations.One issue with using the NPT ensemble for determining the melting point,and particularly for the two-phase simulations, is that the specific temperature ofthe phase boundary proved to be somewhat sensitive to the barostat parameters,and could vary by as much as ±200 K. The barostat sensitivity has been notedpreviously for water simulations107 and was minimized by using larger system sizes.The results given here for the two-phase simulations were obtained employing ananisotropic Parrinello-Rahman barostat103,104 with a relaxation time of 50.0 ps anda compressibility of 3.0×10−5 bar−1 along the three primary axes and 0 otherwise.The slower relaxation time of the barostat dampens volume fluctuations allowing thetwo-phase system to evolve more smoothly into a single phase. The two-phase sim-ulations provide smoother melting point trends, but due to the barostat sensitivity,do not improve the absolute accuracy of the melting point. However, importantly,melting point trends across the model salts were consistent, regardless of the baro-stat settings.39We emphasize that all three temperatures from the hysteresis method (T−,T+, and Tm) as well as Ts and Tl from the two-phase simulations give similar meltingpoint trends, so we would expect the same trends to apply to the thermodynamicmelting temperatures as well.All of the results presented in this thesis are from Gromacs simulations. Iwrote a molecular dynamics (MD) program for this project, but chose to use theGromacs package for the MD simulations for efficiency and convenience. Details ofthe configuration representations and evaluation are given in Appendix A. A briefoverview of how molecular dynamics works is given in Appendix B.2.10 Crystal StructuresExperimentally, some ionic liquids exhibit polymorphism, where the sub-stance can crystallize into multiple crystal structures.84,108 Some of the simple modelsalts considered in this work also exhibit polymorphism, and the melting tempera-tures (at P = 1 bar) of the different crystals can vary by hundreds of degrees.We are interested in which crystal structures are stable near the meltingtransition. When Ts is different for different crystal structures of the same salt,we assume that the crystal structure with the highest Ts is the most stable attemperatures near the melting point. This assumption is reasonable near the meltingpoints because the isotropic liquid is an accessible, equilibrium reference state forboth solids, one of which is more stable than the liquid, the other which is lessstable. Comparing the solids using the liquid as a reference state makes sense attemperatures near the normal melting point. However, this comparison does notallow us to make any statements about which structures are thermodynamically40stable at lower temperatures because the temperature dependence of the chemicalpotential is different for each crystal structure. It is also impossible to discern themore stable structure if the polymorphs have the same Ts with our grid spacing of50 K.The different polymorphs we observe near the melting point are usually dis-tinguished by the degree of orientational order/rotational motion. In the modelswe are considering, we do see salts that likely undergo a s-s transition upon heat-ing from a highly orientationally ordered solid to an orientationally disordered solidbefore melting, provided the pathway barriers (both potential and kinetic) are notprohibitive.With the exception of the references to the actual NaCl and CsCl salts onPage 1, references to CsCl and NaCl are to the archetypal crystal structure. Forsimplicity, we refer to interpenetrating simple cubic lattices of cations and anionsas “CsCl” solids. Similar definitions (with the appropriate ion lattices) are alsoused for “NaCl” solids, “NiAs” solids, and “wurtzite” solids. The archetypal crystalstructures are for monoatomic ions. When the multi-site cations studied in thisthesis are in the arrangement prescribed by the archetypal crystal structure, thecations may be orientationally ordered or disordered.41Chapter 3How Distributed Charge Reducesthe Melting Points of Model IonicSalts ∗3.1 OverviewIn this chapter, we focus on how cation charge distribution reduces the melt-ing points of simple model salts composed of size-symmetric, spherical, monovalentions. We vary the cation geometries with up to one off-center charge site (2L or 1Cgeometry). The anions in this chapter have a single interaction site (1C geometry).The redistribution of charge away from the cation center can reduce the meltingpoint by up to 50% compared with the charge-centered case. Our results demon-strate how charge distribution can reduce the melting point of molecular ionic solidsto temperatures well below those of typical of inorganic salts. The model salts canbe viewed as an interesting step away from inorganic molten salts, toward molecularionic liquids.∗A version of this chapter has been published. E. K. Lindenberg and G. N. Patey, “Howdistributed charge reduces the melting points of model ionic salts,” J. Chem. Phys. 140, 104504(2014).423.2 Melting Point ResultsThe 1C - 1C SaltIt is convenient to begin with the 1C - 1C salt, which serves as a referencesystem from which to measure the influence of charge distribution. The 1C - 1Csalt persists as a supercooled liquid at 1050 K (T−) and solidifies at 1000 K. Thefreezing transformation is signalled by a sudden drop in the average configurational(potential) energy, U¯ , shown in Fig. 3.1. The crystalline solid persists at 1400 K(T+). The melting of the superheated solid is signalled by a sudden increase in theaverage configurational energy between 1400 and 1450 K.43800 1000 1200 1400 1600Temperature (K)-520-500-480-460-440-420U (kJ/mol)CoolingHeatingTwo Phase SimulationTm (Eq. 2.5)T-T+SolidLiquid1C - 1CLc = 0.0TlTsFigure 3.1: The average configurational energies from cooling and heating of the1C - 1C salt. The supercooled liquid freezes between 1050 K and 1000 K, andthe superheated solid melts between 1400 K and 1450 K. The estimated meltingtemperature is 1238 K using Eq. (2.5), and between 1250 K and 1300 K usingtwo-phase simulations.44The melting temperature of the 1C - 1C salt is 1238 K, as estimated byEq. (2.5), and the two-phase simulation result agrees with this estimate. At 1250K (Ts), the two-phase system evolves into a solid and at 1300 K (Tl), the systemevolves into a liquid. The stable crystal structure for the size-symmetric 1C - 1Csalt is a CsCl crystal, which is predicted by both the phase diagram of the RPM109and the radius ratio rules. In the RPM, the fluid and CsCl crystal coexist atT ∗ = kTσ/q2 = 0.0227 and P ∗ = Pσ4/q2 = 0, where  is the dielectric constant ofthe medium, σ is the hard sphere diameter, q is the ion charge divided by√4pi0,and k is the Boltzmann constant.109 For our model parameters (taking σ = 0.50nm and  = 1.0), the RPM solid-liquid (s-l) transition temperature would be 759 Kat P = 0 bar. For a pure LJ system with our model parameters, the s-l transitionwould be near 310 K.110 The 1C - 1C salt has a higher melting temperature thanboth the corresponding RPM and pure LJ systems, as one would expect.The 2L - 1C SaltsWith the s-l phase behaviour of the 1C - 1C salt established as a referencecase, we turn to the 2L - 1C set of salts. Relevant simulation results for the 1C- 1C salt are included in our discussion of trends across salt models (Table 3.1).Average configurational energy hysteresis loops for a selection of 2L - 1C salts areshown in Fig. 3.2. We discuss the group of salts that show hysteresis and clear phasetransitions, as in Fig. 3.2(a-g), before addressing the salts that lack obvious signs ofa first-order transition, as in Fig. 3.2(h).Nearly seventy-five percent of the 2L - 1C salts considered in this chaptershow freezing and melting transitions within hysteresis loops. Specifically, all of thesalts with the 2L33 and 2L50 cations, the 2L67 cations with d ≤ 18, and the 2L10045cations with d ≤ 10 exhibit distinct s-l phase transitions. These salts crystallizespontaneously as CsCl solids with the cation and anion centers occupying the Cs+and Cl− sites, respectively. The off-center charges on the cations do not occupylattice sites. For the salts with relatively small charge arms, the off-center chargesites are essentially in free rotation around the cation center, and the crystal isorientationally disordered.111,112 For salts with larger charge arms, the rotationalmotion is hindered. We return to this point in the discussion of reorientationaldynamics below.The melting point trends are easily identified by comparing T+ or T− inFig. 3.2(a-g). The s-l transition shifts to lower temperatures within each columnof Fig. 3.2 (constant d) and each row (constant fc), indicating that both modelparameters influence the melting temperature.In Fig. 3.2, the vertical distance between the two phases at the meltingtemperature corresponds approximately to ∆fusU¯ , which decreases significantly asthe charge arm increases. The average configurational energies of the solids areconsistently near -525 kJ/mol at 450 K. It is the average configurational energies ofthe liquids that decrease.46-550-500-450U (kJ/mol)-550-500-450U (kJ/mol)-550-500-450U (kJ/mol)600 800 1000 1200 1400Temperature (K)-550-500-450U (kJ/mol)600 800 1000 1200 1400Temperature (K)(c)   2L50-08 - 1C(e)   2L67-08 - 1C(g)  2L100-08 - 1C(d)   2L50-16 - 1C(f)   2L67-16 - 1C(h)  2L100-16 - 1C(a)   2L33-08 - 1C (b)   2L33-16 - 1CLc = 0.11 Lc = 0.21Lc = 0.16 Lc = 0.32Lc = 0.21 Lc = 0.43Lc = 0.32 Lc = 0.64Figure 3.2: The hysteresis loops of the 2L - 1C salts with d = 08 and d = 16show that increasing the charge arm decreases the melting point. Increasing eithermodel parameter shifts the s-l transition to lower temperatures. The 2L100-16 -1C salt does not show any hysteresis in panel (h). This salt becomes trapped in aglassy state with a glass transition temperature, as defined in the text, of 400 K.The legend is the same as Fig. 3.1.47Table 3.1 contains numerical results for the model ILs that crystallize sponta-neously. We report the melting temperature calculated from Eq. (2.5), Tm, roundedto the nearest simulation temperature. The Tm values decrease with increasingcharge arm, and range from 1250 K to 700 K. The temperature-dependent proper-ties given in the last six columns are reported at Tm.The average reduced densities ρ¯∗ for the solid and liquid states are calculatedas Nσ3/V¯ , where N is the number of ions, σ is the LJ length parameter, and V¯is the average simulation cell volume in nm3. The reduced average densities of thesolids do not show much variation from the 1C - 1C value of ρ¯∗s = 1.03 as the chargearm increases and Tm decreases. The liquid phase, however, is more sensitive tochanges in charge arm and temperature. The average reduced density of the liquidρ¯∗l increases from 0.78 (1C - 1C at 1250 K) to 0.99 (2L67-18 - 1C at 700 K).The sixth column in Table 3.1 gives U¯ IPs /U IPiso, which is the ratio of the aver-age configurational energy per ion pair in the solid at Tm to the minimum potentialenergy of an isolated ion pair. The U¯ IPs /U IPiso values are larger than 1.0 for all ofthe salts in Table 3.1, which is expected because all of the salts crystallize spon-taneously. As the charge arm increases, U¯ IPs /U IPiso decreases, which indicates that alarger fraction of the crystal energy comes from ion pairing rather than from thepacking arrangement.The values of ∆fusH¯ are calculated as H¯l - H¯s using the single phase simula-tions from the hysteresis cycles. The enthalpy values of the solids calculated fromthe hysteresis cycles will be artificially high due to the random orientation of thecrystal within the simulation cell and crystal stacking faults. The magnitude ofthe increase varies from salt to salt. We estimate that the potential energies areincreased by up to ∼5 kJ/mol compared to a properly oriented crystal. Artificially48high values of H¯s would lead to underestimations of ∆fusH¯. However, as with themelting temperature estimates, it is the trends in ∆fusH¯ across a series of salts thatare instructive, not the specific values themselves.Table 3.1: Properties of the crystallizing 2L - 1C salts near their melting points.The charge arm Lc is included for easy reference. The temperature Tm is cal-culated from Eq. (2.5) and rounded to the nearest simulation temperature. Thequantities in the last six columns are reported at Tm. The average reduced num-ber densities ρ¯∗ of the liquids show a more pronounced increase than the solidsas the charge arm increases. The ratio of average configurational energies per ionpair in the solid to the potential energy of an isolated ion pair, U¯ IPs /U IPiso, andenthalpies of fusion, ∆fusH¯ (kJ/mol ion pairs), both decrease as the charge armincreases. The relative contributions of the Lennard-Jones (LJ) and electrostatic(ES) energies to ∆fusH¯ are given in the last two columns as percentages.Cation Lc Tm ρ¯∗(s) ρ¯∗(l) U¯ IPs /U IPiso ∆fusH¯ ∆fusH¯— (K) — — — (kJ/mol) %LJ %ES1C 0.00 1250 1.03 0.78 1.70 38 40 602L33-06 0.08 1250 1.03 0.78 1.62 37 40 602L33-08 0.11 1250 1.03 0.79 1.58 36 40 592L33-10 0.13 1200 1.00 0.80 1.55 31 41 582L33-12 0.16 1200 0.96 0.81 1.50 25 41 592L33-14 0.19 1150 1.01 0.83 1.49 28 44 562L33-16 0.21 1150 1.05 0.83 1.46 32 46 542L33-18 0.24 1100 0.99 0.84 1.40 21 48 522L33-20 0.27 1050 1.03 0.86 1.37 23 52 472L33-22 0.29 1050 1.05 0.86 1.33 25 57 422L33-24 0.32 900 1.05 0.90 1.30 19 65 352L50-06 0.12 1150 0.99 0.82 1.58 28 42 582L50-08 0.16 1200 1.04 0.80 1.59 35 42 582L50-10 0.20 1150 1.01 0.83 1.50 29 44 562L50-12 0.24 1150 1.04 0.83 1.46 31 46 54Continued on next page. . .49Table 3.1: (Continued)Cation Lc Tm ρ¯∗(s) ρ¯∗(l) U¯ IPs /U IPiso ∆fusH¯ ∆fusH¯— (K) — — — (kJ/mol) %LJ %ES2L50-14 0.28 1150 1.08 0.83 1.43 34 49 512L50-16 0.32 1000 1.04 0.88 1.37 21 55 452L50-18 0.36 900 1.03 0.92 1.33 16 62 382L50-20 0.40 850 1.06 0.93 1.29 16 65 352L50-22 0.44 800 1.05 0.93 1.26 15 64 362L50-24 0.48 750 1.05 0.93 1.22 15 64 362L67-06 0.16 1250 1.07 0.79 1.56 41 42 582L67-08 0.21 1200 1.03 0.81 1.50 33 43 572L67-10 0.27 1150 1.04 0.83 1.45 30 46 542L67-12 0.32 1000 1.04 0.89 1.41 22 53 472L67-14 0.37 950 1.05 0.90 1.36 19 58 412L67-16 0.43 800 1.05 0.96 1.32 12 68 312L67-18 0.48 700 1.11 0.99 1.29 13 70 302L67-20 0.53 700 1.10 0.97 1.24 15 66 342L100-06 0.24 1100 1.03 0.85 1.49 27 45 552L100-08 0.32 1050 1.03 0.87 1.43 23 50 502L100-10 0.40 1000 1.11 0.89 1.39 26 58 42From the values in Table 3.1, we see that increasing the charge arm reducesboth the melting point and ∆fusH¯. The breakdown of the LJ and electrostatic con-tributions to ∆fusH¯ are given in the last two columns of Table 3.1. For the 1C - 1Csalt, the LJ and electrostatic contributions to ∆fusH¯ are 40% and 60%, respectively.The stabilization provided by the LJ energy is one of the thermodynamic drivingforces for crystallization, even though it accounts for only ∼6% of the average po-50tential energy of the 1C - 1C salt in the CsCl solid. As the charge arm increases,the relative contribution of the LJ interactions to ∆fusH¯ increases, although ∆fusH¯decreases in magnitude. The 2L67-18 - 1C salt has a ∆fusH¯ = 13 kJ/mol, with 70%coming from LJ interactions and 30% from electrostatics. The marked decrease inthe electrostatic contribution to ∆fusH¯ suggests that the configurations in the liquidare comparable, at least in terms of electrostatic energies, to those in the CsCl solid.The LJ contribution to ∆fusH¯ also decreases in magnitude as the charge armincreases, but it contributes an increasing percentage of the ∆fusH¯ value. The 2L67-16 - 1C salt has the lowest ∆fusH¯ of the spontaneously crystallizing salts with a valueof 12 kJ/mol, with 8 kJ/mol from LJ interactions and 4 kJ/mol from electrostaticinteractions. For this salt, the ion pair energy is about 75% of the average configu-rational energy per ion pair in the solid. The majority of the configurational energydifference between the CsCl solid and the liquid for the salts with large charge armscan be attributed to the difference in LJ interactions.The LJ potential energies in the CsCl solids increase (become less negative)with increasing d and, in the case of the 2L50-24 cation, become positive. The strongelectrostatic interaction causes the ions to “penetrate” the repulsive LJ sphere andcreate bound ion pairs that exist in both the solid and liquid phases. The ion pairsare bound tightly and each ion sphere is deformed (flattened) on one side.The values of both U¯ IPs /U IPiso and ∆fusH¯ decrease as the charge arm increases.Both quantities measure the relative stability of the solid, the former with respectto a collection of isolated ion pairs, and the latter with respect to the liquid phase.As the charge arm increases, the salts clearly have less thermodynamic incentive tocrystallize.The molar entropy of fusion (∆fusS = ∆fusH¯/Tm) can be estimated from51the ∆fusH¯ and Tm values given in Table 3.1. For the systems included in Table 3.1,∆fusS also decreases with increasing charge arm, indicating that there is less entropicgain associated the s-l transition. In isolation, the reduction of ∆fusS would usuallyincrease melting temperatures, however, the fractional change in the entropy offusion with respect to the 1C - 1C system is less than that for ∆fusH¯, consistentwith the reduced melting points observed.We have discussed the enthalpies of fusion at length in this section, whichdepend on both the charge arm and the melting temperature. We will revisit theenthalpic properties of the salts again in Chapter 3.4 when we examine the enthalpiesat set temperatures in both the solid and liquid phases.The “Glassy” 2L - 1C SaltsWe now discuss the 2L - 1C salts that do not crystallize spontaneously whencooled, at least not on time scales currently convenient for simulation. Instead,they become trapped in “glassy” states. Because glass formation can depend on thecooling rate, all of the salts that formed glasses when cooled in 50 K increments (tenin total), were again cooled in 25 K increments. Only one of the ten salts (2L67-20- 1C) crystallized at the slower cooling rate, and can be regarded as lying at theboundary between salts that spontaneously crystallize, and those that tend to formglasses. The nine salts that did not spontaneous crystallization include the cations2L67-22, 2L67-24, and 2L100 with d ≥ 12. The glass transition temperatures Tg(defined in Chapter 2.7) for these salts occur between 350 and 450 K, which is nearthe conventional threshold value of 373 K that is used to “define” ionic liquids, asopposed to molten salts.52In order to explore possible crystal structures for the glassy salts, and to ob-tain T+ values for the melting point estimates, we heated prepared crystal structuresfrom T = 50 K. Upon heating, two of the crystals of the 2L100-16 - 1C salt rear-ranged into a highly orientationally ordered, lower enthalpy structure. The spacegroup of the new crystal is 111 (P4¯2m), as identified using the FINDSYM113 program.We refer to this crystal structure as 111n, following the labelling convention used byFilion and Dijkstra114 of giving the space group and occupied Wyckoff position(s).The 111n crystal structure, shown in Fig. 3.3, has a tetragonal unit cellwith two equal axes and one longer axis. There are staggered columns of ion pairsparallel to the long axis of the unit cell. The orientations of the cations are thesame within each column, and alternate from column to column. Perpendicular tothe columns, the crystal structure has alternating layers of electrostatic- and LJ-dominated interactions. The view in the [1¯1¯0] direction, shown in Fig. 3.3(b), depictsthe staggering of adjacent columns and highlights the electrostatic-dominated layers(red and white) and the LJ-dominated layers (blue). The 111n crystal structureresembles the layered crystal structures used in intercalation studies and potentialsuperconductors.115,116 The crystal plane shown in Fig. 3.3(b) bears a resemblanceto iron oxychloride (FeOCl, space group 59, Pmmn), where the Fe, Cl and O positionsare occupied by the cation off-center charge site, cation center, and anion center,respectively.117 The differences between the two crystal structures are visible in theother crystal planes.The enthalpy of the 111n crystal is lower than that of the CsCl crystal for allglassy salts. Upon heating, the 111n structure did not rearrange to another crystalstructure, or show any reorientational relaxation prior to melting, for any of theglassy salts. The orientationally disordered CsCl (OD-CsCl) and 111n structures53are two potential crystal arrangements for our glassy salts near their melting tem-peratures. Our simulations do not indicate which is the more stable crystal. TheOD-CsCl structure is perhaps partially stabilized by the reorientational entropy ofthe cation, whereas the 111n structure is stabilized by a lower configurational energy.(a) 111n unit cell (b) [110] view(d) [001] view(c) [010] viewzx yFigure 3.3: The 111n crystal structure shown for the 2L100-16 - 1C salt. Thecrystal space group is 111 (P4¯2m). The cation centers, cation off-center charge sites,and anion centers are shown in blue, white, and red, respectively. The ions arenot shown to scale; the off-center charge sites on the cations are embedded withinthe cation volume. In (d), the top layer of anions is not shown to highlight thealternating orientations of the cations.54The 2L67-20 - 1C salt lies between the salts that crystallize in a CsCl struc-ture and the salts that vitrify. The salt crystallized during the heating part of thehysteresis loop, shown in Fig. 3.4. We cooled the liquid at a slower rate (25 Kincrements) and the 2L67-20 - 1C salt does crystallize into a CsCl structure. The2L67-20 - 1C salt was the only glassy salt that crystallized at the slower coolingrate. For the CsCl solid, we take T− as 550 K, rather than using the glass transitiontemperature of 400 K. For the 111n solid, we use the glass transition temperatureas T−. Both superheated crystals melt at 900 K.400 600 800 1000Temperature (K)-600-580-560-540-520U (kJ/mol)CoolingHeatingTm (Eq. 2.5)Heat 111n2L67-20 - 1CLc = 0.53TgFigure 3.4: Hysteresis loops for the 2L67-20 - 1C salt. The salt freezes into aglassy state when cooled in 50 K increments, but crystallizes into a CsCl structurewhen it is cooled at a slower rate (25 K increments). The green squares mark theheating of the 2L67-20 - 1C salt in the 111n crystal structure. The 111n structure ishighly orientationally ordered and the cations do not show significant reorientationalrelaxation until the salt melts.55The configurational energies as a function of temperature for the 2L100-16- 1C salt are shown in Fig. 3.5, as a representative example of the glassy salts.The cooling and heating traces of the liquid/glassy systems nearly coincide. Whenheating the prepared crystal structures, the solids rearrange into one of two crystalstructures before melting. The NiAs crystal (not shown) rearranges into the OD-CsCl structure when heated. The OD-CsCl crystal persists between 650 K and 800K, and liquefies at 850 K. The wurtzite and NaCl structures rearrange into the 111nstructure when heated. The 111n crystal melts at 1050 K.400 600 800 1000Temperature (K)-600-580-560-540U (kJ/mol)CoolingHeatingHeat CsClHeat NaClHeat wurtziteHeat 111n2L100-16 - 1CLc = 0.64TgliquidOD-CsCl111 nFigure 3.5: The average configurational energy hysteresis loop for the 2L100-16 -1C salt. When cooling from the liquid state, this salt becomes trapped in a glassystate. The crystal structures heated from T = 50 K rearrange into an orientationallyordered 111n structure or into an orientationally disordered CsCl crystal beforemelting.56(a) Liquid(b) OD-CsCl (c) 111nFigure 3.6: Configurational snapshots of the 2L100-16 - 1C salt in three differentphases at T = 700 K. From the two-phase simulations, Ts for the 111n structureof this salt is 800 K. The cation centers, cation off-center charge sites, and anioncenters are shown in blue, white, and red, respectively. The ions are not shownto scale; the off-center charge sites on the cations are embedded within the cationvolume.57Configurational snapshots of the 2L100-16 - 1C salt at T = 700 K are shownin Fig. 3.6. The three phases included (T = 700 K) are the supercooled liquid, theOD-CsCl solid, and the 111n solid. The two-phase simulations with the 111n crystalgive Ts = 800 K for this salt. The simulation cells are not shown to scale. Thenumber of particles in each simulation varies (liquid N = 432, CsCl N = 686, 111nN = 576).Table 3.2 contains the melting temperatures [from Eq. (2.5)] for both theCsCl and 111n structures, and the reduced densities for the 2L - 1C glassy salts.As in Table 3.1, Tm is rounded to the nearest simulation temperature. The reduceddensities are reported at the Tm of the 111n structure. The temperature range inthe table is only 250 K, so the variations values for the CsCl crystals continue thetrends from Table 3.1, with Tm generally decreasing with increasing Lc. For the2L100 - 1C series, the Tm values of the 111n structure peak between d = 14 and 20and are higher than the corresponding CsCl estimates.In two-phase simulations including the 111n crystal and liquid, the crystalwas placed with the (110) plane facing the liquid. The crystal phase was able togrow and, at certain temperatures, the entire liquid phase crystallized. The lowestvalue of Ts obtained from our set of model salts was 600 K for 2L100-20 - 1C. The2L100 - 1C salts with d = 22 and 24 were exceptions in that two-phase simulationsdid not result in freezing, even at 400 K. It is possible that these two salts have adifferent underlying crystal structure. In any case, for discussion purposes we takeTs to be 400 K for these salts.58Table 3.2: Properties of the 2L - 1C salts that vitrified. The values of Tm are givenfor both CsCl and 111n solids. The reduced average densities are reported at Tm ofthe 111n solid.Cation Lc Tm Tm ρ¯∗CsCl 111n (111n) (l)(K) (K)2L67-20 0.53 700 650 1.084 0.9912L67-22 0.59 650 650 1.079 0.9672L67-24 0.64 650 650 1.064 0.9462L100-12 0.48 900 700 1.073 1.0072L100-14 0.56 750 750 1.093 0.9862L100-16 0.64 650 750 1.114 0.9772L100-18 0.72 600 750 1.117 0.9522L100-20 0.80 650 750 1.097 0.9122L100-22 0.88 600 650 1.088 0.9332L100-24 0.96 — 500 1.108 0.9983.3 Melting Point TrendsOur simulations give five transition-related temperatures for each salt: T−,T+, Tm from Eq. (2.5), and Ts and Tl as the boundaries of the two-phase simulations;all are presented concisely in Fig. 3.7. All five transition-related temperatures showa general decrease as the cation charge is moved further off-center. Since the qual-itative behaviour of all five temperatures is the same, we would expect the generaltrends to apply as well to the true thermodynamic melting temperatures.598 12 16 20 24 8 12 16 20 240400600800100012001400Temperature (K)T- to T+Ts  to Tl Tm (Eq. 2.5)Tg to T+ (111n)T+ (CsCl)400600800100012001400Temperature (K)1C 2L33 2L502L67 2L100d value (0.01 nm)1C***Figure 3.7: Melting temperature summary for the 2L - 1C salts. All five of thetransition-related temperatures produce the same trends. The melting point de-creases with increasing charge arm with the exception of the 2L100 cations withd = 12 − 16. The increase in melting points coincides with the change in crystalstructure from CsCl to 111n. Tabulated values are given in Appendix D in TablesD.1 and D.2.60As discussed in Chapter 2, the two-phase simulations were performed torefine the melting temperature estimates, and here we use Ts to discuss particulartrends. In order to highlight the differences in melting points, we scale the Ts valueof each salt by Ts of the reference 1C - 1C salt (1250 K).The plots in Fig. 3.8 show how the melting point changes when one ionmodel parameter, fc or d, is varied holding the other fixed. In Fig. 3.8, we haveplotted only the melting points for the salts that spontaneously crystallize into aCsCl structure. The melting points generally decrease as fc or d increases, but thedecrease is not uniform.0.2 0.4 0.6 0.8Off-center charge (e) / Ts(1C - 1C)060810121416182022240.2 0.4 0.6 0.8d/Rion2L33 - 1C2L50 - 1C2L67 - 1C2L100 - 1C(a) (b)Figure 3.8: The two parameters of the 2L model ions, the amount of charge off-center, fc, and the displacement distance from the ion center, d, both influence themelting temperatures of the 2L - 1C salts.61The main results of this chapter are summarized in Fig. 3.9, where the ratioof Ts/Ts(1C - 1C) is plotted versus the cation charge arm. The influence the chargearm has on melting point is clear; moving some of the unit charge away from thecenter in an RPM-like salt generally decreases the melting temperature. For thesalts near Lc = 0.50, where both CsCl and 111n structures are viable possibilities,we plot the highest Ts value observed for each model salt in Fig. 3.9. Fig. 3.9 rep-resents the minimum reduction in melting point for our model salts.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1Normalized Charge Arm, Lc0. / Ts(1C - 1C)2L33 - 1C2L50 - 1C2L67 - 1C2L100 - 1COD-CsCl 111nFigure 3.9: The plot of scaled melting temperatures against the normalized chargearm reveals that melting points of CsCl crystals decrease until the value of Lc reaches0.50. Above Lc = 0.50, the 111n has a lower enthalpy than the CsCl structure.The 111n melting temperature peaks between Lc = 0.64 and 0.72. All of saltswith Lc ≥ 0.50 vitrify; none of these salts spontaneously crystallize into the 111nstructure from the liquid phase on simulation time scales.62The salts with Lc ≤ 0.50 crystallize in a CsCl structure. For equivalentcharge arm values below Lc = 0.33, the 2L33 - 1C salts melt at temperatures equalto, or lower than, the other model salts. The 2L33 - 1C salts have a maximum Lcvalue of 0.33. For 0.33 < Lc ≤ 0.50, the 2L50 - 1C salts are the lowest melting CsClsalts, until the maximum Lc = 0.50 is reached. For a given Lc a model cation witha larger d achieves a lower melting point than one with a larger amount of chargeoff center. This feature may be useful for the rational design of low-melting ILs.At charge arms between 0.50 and 0.64, we see different trends for the 2L67- 1C salts and the 2L100 - 1C salts. In this region, the highest melting crystaltransitions from CsCl to 111n. The two model salts in the crossover region are the2L67-20 - 1C and 2L100-12 - 1C salts. The melting points of the CsCl solid werehigher than those of 111n for both of these salts. As the charge arm increases further,the melting points of the 2L67 - 1C salts continue to decline, despite the change incrystal structure, whereas the melting points of the 2L100 - 1C salts increase. Themelting point increase in the 2L100 - 1C salt series occurs between d = 14 and 16.At d = 14, the ion geometry and the 111n structure are compatible, but the optimalion geometry for the 111n structure occurs at larger charge arms (d = 16 to 18). Thepeak in the melting points lies near 800 K for the 2L100 cations with 0.64 ≤ Lc ≤0.72. Once the charge arm is above 0.72, the melting points drop again.We have shown that varying the charge distribution in size-symmetric, spher-ical, monovalent salts can lower the melting temperature by more than 50%. Uponcooling, the model salts considered here form orientationally disordered crystalsor vitrify, two solidification behaviours that have been observed in real IL sys-tems.48,68,75 In the next section, we analyze thermodynamic, structural, and dynamicproperties in both phases to better understand why the melting point decreases with63increasing charge arm.3.4 2L - 1C Salt PropertiesThe thermodynamic melting point occurs when the chemical potentials ofthe solid and the liquid are equal at the same temperature and pressure. Typically,a solid is stabilized by enthalpy and a liquid by entropy. The average enthalpy iseasily accessible in MD simulations, and it is the first property we examine. Thecharge arm alters both the symmetry and strength of the electrostatic interactions,and impacts the average enthalpy in both the solid and liquid phases.The average enthalpies at T = 1100 K, shown in Fig. 3.10, offers insightinto one of the main reasons why melting temperatures decrease with increasingcharge arm. The horizontal black lines give the average enthalpy of the 1C - 1Csalt in the solid phase (solid line) and supercooled liquid phase (dashed line) atT = 1100 K. The plot shows how the charge arm affects the average enthalpy ofeach phase as d increases. The ions with d = 6 are on the left side of Fig. 3.10, andthe enthalpy values deviate slightly from those of the 1C - 1C reference salt. Theenthalpies of the CsCl solids and the supercooled liquids at d = 6 are slightly higher(less negative) than the enthalpies of the reference 1C - 1C salt, but the enthalpicgap between phases remains near 40 kJ/mol. As d increases, the average enthalpydecreases for both the supercooled liquid and CsCl solid phases, primarily due tothe increasing attractive strength of the electrostatic interaction from the off-centercharge site. The rate of decrease is faster for the supercooled liquid, which indicatesthat the liquid phase benefits more from the charge redistribution than the CsClsolid. For each set of salts, there is a d value where the average enthalpy of the64solid increases. This increase marks where the CsCl solid becomes unstable withrespect to the liquid across salt models. To the right of the increase, the salts aremolten and the average enthalpies follow the trace of the liquid. At large d values,the average enthalpies of the liquids are significantly lower than those of the CsClsolids. The increase in charge arm stabilizes the liquids more than the CsCl solids,and the melting points shift to lower temperatures.656 8 10 12 14 16 18 20 22 24d (0.01 nm)-500-475-450-425H (kJ/mol)1C2L332L502L672L100liquidssupercooled liquidsCsCl solidsT = 1100 KFigure 3.10: The average enthalpies as a function of charge displacement at T =1100 K reveal one of the main reasons the melting points decrease. Moving fromleft to right, the enthalpies of the supercooled liquids (marked with dot-dash lines)decrease faster than the enthalpies of the solids (marked with solid lines) as thecharge arm increases. The single enthalpic increase in each series marks the end ofthe superheated CsCl phase, and the approximate phase transition at 1100 K. Onthe right, the 2L50, 2L67, and 2L100 series salts are liquids (marked with dashedlines) that have significantly lower enthalpies than the reference 1C - 1C CsCl solid(solid black line).66A similar plot of the enthalpies at T = 600 K, shown in Fig. 3.11, showsthat the 111n structure has a higher enthalpy compared to the CsCl solid for the2L100-10 - 1C salt, and a lower enthalpy for the 2L100-12 - 1C salt. The crossoverbetween crystal types coincides with the boundary between salts that spontaneouslycrystallize and those that do not. For 2L67 cations with d ≥ 20, the 111n structureshave a lower enthalpy than the CsCl structures. The 2L67-20 - 1C salt is on theboundary between the solidification behaviours but does crystallizes spontaneouslyinto an OD-CsCl structure. The enthalpy traces across salt models still show thetransition when the liquid becomes more stable than the solid but, in contrast toFig. 3.10, the transition does not involve a marked increase in enthalpy.In Fig. 3.10 and Fig. 3.11, the average enthalpies for the solids were calculatedusing prepared crystals where the orientation of the crystal axes and simulation cellaxes are aligned. The axes alignment gives smoother trends for discussion.676 8 10 12 14 16 18 20 22 24d (0.01 nm)-750-700-650-600-550-500H (kJ/mol)   1C   2L67   2L100             `   CsCl   Liquid   111nT = 600 KLine style   PhaseFigure 3.11: Average enthalpies as a function of charge displacement at T = 600 K.The line styles denote the phase and the symbol shapes/colours denote the cationtype. The solid black line is the enthalpy of the 1C - 1C CsCl solid. The CsCl solidsfor the 2L - 1C salts are shown as solid lines. For 2L100 cations with d ≥ 12 and2L67 cations with d ≥ 18 the (supercooled) liquid enthalpies are connected withdashed lines. The dot-dash-dot lines connect the enthalpies of the salts in the 111nstructure.68For the remainder of this section, we highlight the changes in structural anddynamical properties by focusing on a subset of the 2L67 - 1C salts. The 2L67 cationshave 2/3 of the unit charge located on the off-center site. The 2L67 cations withd = 06, 10, 14, 18, and 22 capture the typical behaviour with respect to changesin charge arm and temperature. It is worth remarking that in real ILs, at lowtemperatures cations with shorter side chains show qualitatively similar behaviourto cations with longer side chains at higher temperatures.12 We see qualitativelysimilar behaviours between salts with small cation charge arms at low temperaturesand salts with large cation charge arms at higher temperatures. Both the alkyl chainlength and charge arm represent at some level the degree of displacement betweenthe center of charge and the center of mass.We discuss properties of the salt subset at 500 K and 1100 K. The salts andtheir phases at T = 500 K and 1100 K are given in Table 3.3.Table 3.3: The phase behaviour of the 2L67 - 1C salt subset at T = 500 K and1100 K listed for reference. The dynamic properties of these salts are discussed inthe rest of this section.Temperature 2L67-d - 1C Salt Subset(K) d = 06 10 14 18 22500 CsCl CsCl CsCl CsCl slow liquid or 111n1100 CsCl CsCl liquid liquid liquidThe radial distribution functions (rdfs) at T = 500 K are shown in Fig. 3.12.In panels (a-d), the salts are CsCl crystals. In panel (e), the 2L67-22 - 1C salt isa slow liquid nearing its glass transition temperature and shows little long-rangeorder. In panel (f), the 2L67-22 - 1C salt is shown again, but in the 111n crystalstructure. Across the set of rdfs, the changes in structure due to increases in thecharge arm are most pronounced in the CC-AM rdf (light blue). In panels (a) and69(b), the first peak in the CC-AM rdf is broad and centered around the first CM-AMpeak. The CC sites on the cations are essentially in free rotation around the CMsites. As the charge arm increases, a pronounced shoulder appears on the left ofthe first CC-AM peak in panel (c), and it develops into a separate peak in panels(d-f). The splitting of the first CC-AM peak indicates that the CC sites on thecations prefer to be directed towards particular anions rather than rotating freelyover all orientations. In panel (e), the combined evidence of the position of thefirst peak in CC-AM, and the similarities in the CM-CM (black) and AM-AM (red)rdfs, suggests that the salts are developing ion pairs. In the 111n structure shownin panel (f), the short CC-AM separation signals ion pairing, but the differences inCM-CM (black) and AM-AM (red) rdfs indicate a more complicated arrangementthan simple ion pairs.At 1100 K, the solids in Fig. 3.13(a-b) are still CsCl crystals, and the rdfpeaks show the expected thermal broadening and shifts to longer distances. Theliquids in Fig. 3.13(c-e) have the left shoulder (peak) at shorter separations in theCC-AM rdfs. The 2L67-18 - 1C and 2L67-22 - 1C salts show signs of ion pairingas discussed above. The liquids do not show significant ordering beyond 1 nm. Inthe liquid phase, the ion pairing reduces the longer range charge ordering that canoccur in inorganic molten salts.70246g(r)246g(r)246g(r)246g(r)CM - CMCM - AMAM - AMCC - AM246g(r)0 0.5 1 1.5r (nm)246g(r)(e)           2L67-22 - 1CLc = 0.16solid(a)   2L67-06 - 1C(b)   2L67-10 - 1CLc = 0.27solid(c)   2L67-14 - 1CLc = 0.37solid(d)   2L67-18 - 1CLc = 0.48solidLc = 0.59liquid(f)           2L67-22 - 1CLc = 0.59solidFigure 3.12: Radial distribution functions for selected 2L67 - 1C salts at T = 500K and P = 1 bar. The CC-AM rdf (light blue, cation off-center site to anion centralsite) shows the largest variation as d increases. In panels (a) and (b), the cation isin nearly free rotation, whereas in panels (e) and (f), the cation is part of an ionpair and cannot rotate freely. The solids in panels (a-d) are OD-CsCl structures,and in panel (f) the salt has a 111n structure.71246g(r)246g(r)246g(r)0 0.5 1 1.5r (nm)246g(r)CM - CMCM - AMAM - AMCC - AM246g(r)(e)           2L67-22 - 1CLc = 0.16solid(a)   2L67-06 - 1C(b)   2L67-10 - 1CLc = 0.27solid(c)   2L67-14 - 1CLc = 0.37liquid(d)   2L67-18 - 1CLc = 0.48liquidLc = 0.59liquidFigure 3.13: Radial distribution functions for selected 2L67 - 1C salts at T = 1100K and P = 1 bar. In panels (a) and (b), the salts are OD-CsCl solids, whereas inpanels (c-e) the salts are liquids. The liquids show little long-range order. As inFig. 3.12, as the value of d increases the CC-AM rdf (light blue) shows the largestvariation in the shape and location of the first peak.72Ion pairing was observed in earlier studies of the dynamical properties ofsimilar model salts in the liquid phase.60,61 Further evidence of ion pairing is foundin the velocity autocorrelation functions (VACFs). Normalized VACFs are definedin the usual way,118–120Cv(t) =〈v(t) · v(0)〉〈v(0) · v(0)〉, (3.1)where v(t) is the velocity at time t.In Fig. 3.14, we show VACFs for three of the 2L67 - 1C salts at T = 1100K; the behavioural trends at 500 K are similar. In Fig. 3.14(a), the 2L67-06 - 1Csalt is a CsCl crystal. The negative regions indicate caging effects, which all threeinteraction sites exhibit. In panel (a), the plots for the ion centers, CM and AM, aredistinct, indicating that the ion centers are moving somewhat independently. TheVACF of the CC site follows the AM site; the CC and AM sites move nearly in phasewith the AM site having a slightly larger amplitude. In panel (b), the VACFs of theCM and AM sites are similar, and the CC VACF has a few oscillations about zero.In panel (c), the CC VACF oscillates about zero and the plots for the ion centers,CM and AM, nearly coincide. For the 2L67-22 - 1C salt, the ions are clearly boundin pairs and are not moving independently. The rapid oscillation in the CC VACFsuggests that the ion pair is undergoing librational/vibrational motion. The VACFsin panel (c) provide dynamical evidence of strong ion pairing in the liquid phase.Further insight into ion-ion correlations and the nature of the solid and liq-uid phases can be obtained by considering reorientational autocorrelation functionsof the cation. For example, if ion pairs are formed we would expect the cationsto exhibit longer reorientational correlation times. Reorientational autocorrelation73-0.500.51VACF0 0.25 0.5 0.75 1time (ps)-0.500.51VACF-0.500.51VACFCMCCAM(a)  2L67-06 - 1C(b)  2L67-14 - 1C(c)  2L67-22 - 1CsolidliquidliquidFigure 3.14: VACFs for selected 2L67 - 1C salts at 1100 K. The CsCl solid in panel(a) shows negative regions, indicating ion caging. In panels (b) and (c), the VACFsfor the ion centers CM and AM are similar and the CC VACF oscillates about zero.In combination, these two characteristics suggest that ion pairs are formed in theliquid phase.74functions are calculated as118–120Cl(t) =1NN∑i=1〈Pl(µi(t) · µi(0)) 〉, (3.2)where Pl is the lth rank Legendre polynomial and µi(t) is the unit vector from CMto CC on cation i at time t. We use the first-rank Legendre polynomial to calculatethe correlation function for the subset of 2L67 - 1C salts at 500 K and 1100 K.Normalized results are shown in Fig. 3.15. At both temperatures, C1(t) decays onpicosecond time scales for cations with d = 6, 10, 14. The decay becomes slower asthe charge arm increases.The negative regions in the reorientational correlation function for the 2L67-06 and 2L67-10 cations at T = 500 K indicate that the cations are behaving essen-tially as free rotors undergoing inertial motion. Both salts are OD-CsCl solids atboth temperatures, and the cations reorient with µ(t) changing direction easily andoften. There is no steric hindrance to prevent the cations from reorienting withinthe CsCl crystal. The “cage” created by the anions is symmetric, and there is littleor no change in the occupied space when the cation reorients. The energetic barri-ers to cation reorientation are low, and the increased entropy from the reorientationincreases the stability of the CsCl crystal.The reorientational correlation function for the 2L67-14 cation does notchange dramatically between T = 500 K and 1100 K, although the salt does gothrough a s-l phase change between these two temperatures. In both phases, aftera brief inertial period, slower decay is observed.The 2L67-18 - 1C salt is a CsCl crystal at 500 K, but the cations are unableto reorient freely. The value of C1(t) reaches 0.25 near 50 ps at T = 500 K, and then750 0.5 1 1.5 2Time (ps)00.51C1(t)00.51C1(t)d = 06d = 10d = 14d = 18d = 22d = 22 (111n)T = 500 KT = 1100 KFigure 3.15: The reorientational correlation functions, C1(t), for selected 2L67 - 1Csalts at T = 500 K and T = 1100 K. At 500 K, the 2L67-06 and 2L67-10 cations(CsCl solids) exhibit negative values indicating inertial rotation. At 1100 K, thesetwo salts remain CsCl solids, and the decay remains on picosecond time scales. Thecorrelation functions for the 2L67-14 cations decay on picosecond time scales at both500 K (CsCl solid) and 1100 K (liquid). At 500 K, the 2L67-18 (CsCl solid) and2L67-22 (slow liquid) cations show little decay on picosecond time scales. At 1100K (liquids), these two salts form ion pairs and exhibit hindered rotational motion.76decays even more gradually. In the CsCl solid at T = 500, the relaxation essentiallyoccurs in “jumps” with µ redirecting from one anion neighbour to another. Themajority of the reorientation in the liquid (at T = 1100 K), occurs within the first3 picoseconds.For the 2L67-22 - 1C salt at T = 500 K, we show reorientational correlationfunctions for a slow liquid, and for the 111n crystal. Both systems show strong, long-lived orientational correlations. The slow liquid decays on nanosecond time scales,whereas all salts in the 111n crystal structure remained orientationally ordered untilmelting. At T = 1100 K, the 2L67-22 cation (liquid) does show relaxation, and thereorientational correlation function decays on the order of tens of picoseconds.In general, the liquids at T = 1100 K show an initial rapid decay, correspond-ing to molecular librations,121 followed by more gradual decay after the first 0.5 ps,which is consistent with a rotational diffusion mechanism. Roy and coworkers alsofound that there were two components to the rotational relaxation in model ILs,the first on subpicosecond time scales and the second on picosecond to nanosecondtime scales.59The reorientational correlation functions corroborate the findings from therdfs and the VACFs. Cations with small charge arms are essentially in free rotation.As the charge arm increases, the ions are more likely to form ion pairs. The ionpair lifetimes increase with increasing charge arm,60 and at large charge arms ionsspend significant time bound in ion pairs.For the salts with large charge arms, the liquid phase does not exhibit spatialor orientational ordering. An orientationally ordered 111n crystal melts to form anorientationally disordered liquid in a single step. Cooling model ILs with large chargearms removes energy from both translational and rotational modes. The gradual77reduction of the rotational diffusion rate on cooling is unlikely to naturally restorethe orientational order present in the 111n crystal. Model salts with large chargearms tend to vitrify likely because the cooperative orientational relaxations requiredto spontaneously form an orientationally ordered molecular crystal are simply muchtoo slow, at least on simulation time scales.3.5 Chapter SummaryWe have shown that, in general, the melting point decreases with increasingcharge arm, and for salts with identical charge arms, the melting point is lowestfor the salts with the largest charge displacement d. These observations may proveuseful for the rational design of low melting salts.We emphasize that the ability of the cation to rotate can strongly influencethe crystal structure, melting point, and behaviour of a particular IL as it is cooled.We would expect the CsCl structure to be the stable crystal for all of the 2L -1C model salts, provided the cations are relatively free to reorient. Cations withsmaller charge arms are relatively free to rotate and crystallize spontaneously intoorientationally disordered CsCl structures. Cations with larger charge arms formlong-lived, directional ion pairs, and the salts do not spontaneously crystallize intoany crystal structure. The model salts tend to become trapped in glassy states onsimulation time scales. For these “glassy” salts, we found an orientationally ordered111n structure that has a lower enthalpy than the CsCl crystal. The low energy111n crystal requires that the cations be orientationally ordered. Upon heating, theorientational order of the cations persists in the 111n structure until the salt melts.We conjecture that the entropy cost of this orientational order creates a nucleation78barrier that is too high for spontaneous crystallization of 111n to be observed onsimulation time scales.The charge arm influences both the enthalpy and entropy of fusion. At leastfor the CsCl crystals, both ∆fusH¯ and ∆fusS decrease with increasing charge arm,the former favouring and the latter opposing melting point reduction. However, thechange in ∆fusH¯ outweighs that of ∆fusS resulting in the reduced melting pointsobserved. Physically, the directional electrostatic forces arising from charge dis-placement lower the average enthalpy of both the liquid and CsCl solid phases, butas a function of charge arm, the decrease is faster in the liquid than in the solid. Thisleads to decreasing values of ∆fusH¯ and lower melting points of the CsCl crystals asthe charge arm increases.79Chapter 4Melting Point Trends and SolidPhase Diffusivities ofSize-Asymmetric Model Salts withDistributed Cation Charge ∗4.1 OverviewIn this chapter, we study the combined effects of charge distribution andsize asymmetry with size-asymmetric 1C - 1C and 2L67 - 1C salts. The ion sizeratios are varied in two ways. In Set A, the distance between cation-anion centersis constrained, σ± = 0.50 nm, whereas in Set B, the anion size is held constant,σ− = 0.50 nm, as discussed in Chapter 2.1. The “size” of the ions is consistentwith common anions in real ionic liquids, in between BF –4 and PF –6 . The resultsin this chapter are divided into four parts. In Chapter 4.2, we briefly discuss thesize-symmetric salts (labelled A100), which act as a reference set for both size-asymmetric sets, and then examine the Set A size-asymmetric salts. In Chapter 4.3,we present the second set of size-asymmetric salts (Set B). In Chapter 4.4, we discussthe solid state structures and diffusion dynamics of the A300 series of salts.∗Parts of this chapter are being prepared for submission as a journal article.804.2 Size Set A (Constrained σ±)The size-symmetric salts were studied in detail in Chapter 3. The chargedisplacement distance on the cation is a key factor in whether or not the cationis able to reorient, which affects the solidification behaviour of the salts. At smalldisplacements, the cation is in essentially free rotation and the salts spontaneouslycrystallize into an orientationally disordered cesium chloride (OD-CsCl) solid. Atlarge displacements, the reorientational motion is hindered such that the ions canform directional ion pairs, and cooling the salts gives “glassy” states rather thancrystals. The underlying crystal structure is orientationally ordered. The coopera-tive orientational relaxations required to spontaneously crystallize are much to slow,at least on current simulation time scales.2 Properties and numerical results for thesize-symmetric salts are included below in Table 4.1.When the ions are size symmetric, there is one uniform length scale for allthree types of ion-ion interactions (cation-cation, cation-anion, and anion-anion).Lifting the condition of size symmetry removes the uniformity and creates threedistinct length scales for each type of ion-ion interaction. In Set A, we changethe ion size ratios while minimizing the changes in cation-anion interactions. Weconstrain all of the salts in Set A to have the same σ± while varying σ+ and σ−.The specific size ratios considered are given in Table 2.1.We first consider the charge-centered 1C - 1C salts in Set A. The seven saltswe consider span size ratios from 1:1 to 3:1. One might expect that the series of1C - 1C salts would follow the predictions from the radius ratio rules for chargedhard spheres and crystallize as CsCl, NaCl, or wurtzite structures. The 1C - 1Csalts from A100 through A233 spontaneously crystallize in CsCl (A100 and A133)81or NaCl (A167, A200, and A233) structures, as predicted. According to the ra-dius ratio rules, the A267 and A300 1C - 1C salts are expected to crystallize in awurtzite structure; instead, these salts form NaCl crystals. We suspect that theLJ contribution to the interaction potential favours the NaCl arrangement with sixneighbouring counterions rather than the wurtzite structure with four. The pre-pared wurtzite crystals of the A267 and A300 1C - 1C salts rearrange into thelower-enthalpy NaCl structure during heating.An example of a hysteresis loop is shown in Fig. 4.1, where the averageconfigurational (potential) energy of the A300 1C - 1C salt is plotted as a functionof temperature. The thermal width of the hysteresis is 450 K, which is typical for thesalts considered here. The average potential energy of the NaCl solid formed uponcooling is higher than the prepared NaCl crystal energy because the spontaneouslyformed crystal axes and simulation cell axes are not aligned and the solid has astacking fault. Both liquid and solid phases persist through the entire two-phasesimulation at T = 1100 K for the A300 1C - 1C salt. The persistence of both phasesindicates that the salt is at or near its melting temperature, rather than clearlyabove or below it. We take the apparent coexistence temperature as Ts and the firsttemperature that gives a pure liquid as Tl.82400 600 800 1000 1200 1400Temperature (K)-500-480-460-440-420U (kJ/mol)CoolingHeatingHeat NaClHeat WurtziteTwo-PhaseSimulation (NaCl)Tm (Eq. 2.5)A300 1C - 1CFigure 4.1: Hysteresis plot of the A300 1C - 1C salt. This salt freezes into a NaClstructure. At a size ratio of 3:1, the expected crystal structure for charged hardspheres is wurtzite. The wurtzite crystal we prepared and heated rearranges intoNaCl before melting. The five transition-related temperatures for this salt are: T− =750 K, T+ = 1200 K, Tm = 1001 K, Ts = 1100, and Tl = 1150 K.83The five transition-related temperatures for the series of 1C - 1C salts in SetA are shown in Fig. 4.2. The melting points of the A100, A200, and A233 1C - 1Csalts are highest, with Ts = 1250 K. These three size ratios coincide with the optimal(or near-optimal) matches between size ratio and packing structure for CsCl andNaCl crystals. The decrease in melting point of the A133 1C - 1C salt marks thedestabilization of the CsCl structure, and the subsequent increase in melting pointof the A167 1C - 1C salt marks the onset of stabilization of the NaCl structure. ForA267 and A300 1C - 1C salts, the melting point decreases which indicates that theNaCl structure is becoming less stable.Increasing the size ratio from 1:1 to 1.33:1 decreases the melting point byabout 100 K. Doubling the size ratio from 1.33:1, to 2.67:1 results in the samemelting point, with a peak at intermediate size ratios (between 2:1 and 2.33:1).Beyond size ratios of 2.33:1, the 1C - 1C salts show decreasing melting point trends.Over size ratios varying from 1:1 to 3:1, Ts decreases by about 150 K which indicatesthat size ratio alone has a marginal effect on the melting temperature.84100 133 167 200 233 267 300Cation Size (% of Anion Size)8001000120014001600Temperature (K)T- to T+Ts to TlTm (Eq. 2.5)1C - 1C SaltsFigure 4.2: Transition temperature summary for the 1C - 1C salts in Set A. Thesalts in Set A isolate the impact of ion size ratio on melting point by keeping thevalue of σ± fixed and varying both σ+ and σ−. The overall reduction in Ts as the ionsize ratios vary from 1:1 to 3:1 is about 150 K. The melting temperature peak forsize-symmetric ions (A100) is for the CsCl structure, while the maximum meltingtemperatures for the NaCl structure occur at size ratios of 2:1 and 2.33:1.85We have established the melting points, and melting point trends, of thecharge-centered salts in Set A. We now use these 1C - 1C salts as starting pointsfor redistributing the cation charge. For each size ratio, we estimate the meltingtemperatures for a series of salts where 2/3 of the charge is moved away from thecation center by a distance d.Example hysteresis loops for four salts are shown in Fig. 4.3. In Fig. 4.3(a),the molten salt vitrifies and the cooling and heating traces nearly coincide. Theconfigurational (potential) energies of the CsCl structure and 111n structure (shownin Fig. 3.3) of this salt increase with increasing temperature at different rates. Thestructure with the lowest potential energies switches from 111n to CsCl between 650K and 700 K. The 111n crystal has a lower melting point than the CsCl crystalbecause the CsCl solid permits cation reorientation while the 111n solid does not.The salts in Fig. 4.3(b) and (c) are typical examples of the hysteresis loops forthe spontaneously crystallizing salts. In both Fig. 4.3(b) and (c), the salts freeze asNaCl solids. In (b) there is a slight rearrangement on heating that produced a singlecrystal, rather than a solid with stacking faults. In Fig. 4.3(d), the salt undergoes thephase transition in stages in both directions. We discuss the premelting structuresand the ion dynamics in more detail in Chapter 4.4.Of the salts in Set A, 38 out of 48 spontaneously crystallize. At size ratiosof 2:1 and larger, all of the salts spontaneously crystallize into NaCl solids. Theten salts that do not spontaneously crystallize were the A100 2L67 - 1C salts withd = 18 and d = 22, and the 2L67 - 1C salts with 14 ≤ d ≤ 26 with size ratios A133and A167.In Chapter 3, the A100 2L67-18 - 1C salt spontaneously crystallizes intoa CsCl structure. Here, with a larger number of ions in the simulation cell, the86300 600 900Temperature (K)-560-540-520-500-480U (kJ/mol)-560-540-520-500-480U (kJ/mol)300 600 900Temperature (K)CoolingHeatingHeat CsClHeat 111nHeat wurtziteTwo-Phase Simulations2L67-14  - 1C 2L67-18  - 1C2L67-18  - 1C2L67-18  - 1C(a) (b)(d)(c)A133 A233A300A267TgpremeltingonsetFigure 4.3: The hysteresis loops of four 2L67 - 1C salts in Set A. In (a), the A1332L67-14 - 1C salt vitrifies upon cooling. The two-phase simulations with each solid,CsCl and 111n, give different melting temperatures, with CsCl melting between 850and 900 K, and 111n melting between 750 and 800 K. This salt is on the borderlinebetween salts with stable CsCl solids and 111n solids. In (b) and (c), the A2332L67-18 - 1C and A267 2L67-18 - 1C salts both spontaneously crystallize as NaClsolids. In (d), the A300 2L67-18 - 1C salt crystallizes and melts in phases.87salt vitrifies. It also vitrifies when cooled in smaller temperature increments of 25K. From heating the set of crystalline solids, we find that the CsCl crystal has thehighest T+. The CsCl structure is likely the solid phase that lies nearest to the liquidfor this salt under ambient pressures, which is consistent with the work describedin Chapter 3. At lower temperatures, the A100 2L67-18 - 1C salt may undergosolid-solid transitions to a structure with more orientational order. At a slightlylarger d, the A100 2L67-22 - 1C salt also vitrifies, which is consistent with theresults in Chapter 3. The underlying crystal structure of this salt is orientationallyordered. One candidate crystal structure is a variation on a CsCl structure withorientational order. Some snapshots and details about the crystal structure are givenin Appendix E. We heated this structure and found T+, but did not do two-phasesimulations with this structure.The A133 and A167 salts that vitrify follow the same pattern as the A100series. As d increases, the melting point decreases and the salts become trappedin glassy states once the reorientational motion of the cations is hindered. Thehysteresis cycle for the A133 2L67-14 - 1C salt is shown in Fig. 4.3(a), and thelow-temperature behaviour of the crystals was discussed above. The A167 2L67-14- 1C salt melts from a NaCl structure. The A133 2L67-18 - 1C and A167 2L67-18 -1C salts both show the CsCl and wurtzite crystals transition to an orientationallydisordered solid before melting. For the A133 2L67-18 - 1C salt, we find that theorientationally ordered 111n solid has the highest melting temperature. However,A167 2L67-18 - 1C salt rearranges into an orientationally disordered NaCl structurebefore melting. For both A133 and A167 size ratios at d ≥ 22, the orientationallyordered solid melts directly, while heating other crystal structures results in lowermelting points and subtle rearrangements before melting.88Table 4.1 gives some of the properties of Set A salts, where the d = 0 entrycorresponds to the 1C - 1C salt. In Table 4.1, Ts gives the highest temperatureof the stable solid from the two-phase simulations for each salt, with the relevantcrystal structure given in the third column. For the polymorphic salts and saltsthat show premelting transitions, we give the highest melting point because we areinterested in the minimum melting point reductions.The temperature-dependent properties of each salt are reported at Ts. Theaverage reduced densities, ρ¯∗, for the solid and liquid states are calculated as (N+σ3++N−σ3−)/V¯ , where N and σ refer to the number of ions and the LJ length parameterwith subscripts referring to the cations (+) or anions (−), V¯ is the average simula-tion cell volume in nm3. The ratio U¯ IPs /U IPiso is a measure of how much stability thecrystal structure gains relative to the minimum potential energy of an isolated ionpair, with higher numbers indicating that there is more thermodynamic incentive forthe salt to crystallize. The enthalpy of fusion, ∆fusH¯, at Ts is estimated as H¯l− H¯s,where the average enthalpies of each phase at Ts are calculated from single-phasesimulations.The trends within most size ratios follow the trends of the size-symmetricsalts. The average reduced densities of the liquids generally increase while Ts,U¯ IPs /U IPiso, and ∆fusH¯ generally decrease with increasing cation charge displacement.For the 1C - 1C salts in Set A, U¯ IPs /U IPiso decreases monotonically as the size ratioincreases, from 1.72 to 1.62. In Set A, the heats of fusion range from 44 kJ/mol ofion pairs (A100 1C - 1C) to 6 kJ/mol of ion pairs (A300 2L67-26 - 1C). The s-l tran-sition occurs in stages for the A300 2L67-26 - 1C salt; the low ∆fusH¯ is a reflection ofthe increased solid enthalpy after the premelting transitions at lower temperatures.Another indication that the solids of the A267 and A300 salts undergo premelting89transitions comes from a notable decrease in the average reduced densities of thesolid as d increases. The premelting transitions are discussed in Chapter 4.4 below.Table 4.1: Properties of the Set A salts. The crystal structures marked with anasterisk likely have a s-s transition to a more orientationally ordered crystal struc-ture at lower temperatures. The var-CsCl structure is shown in Appendix E.2.The crystal structures listed as fcc(+) and fcc(PC) melt from a cation fcc latticeand fcc plastic crystal phase, respectively.Size d Crystal Ts ρ¯∗(s) ρ¯∗(l) U¯ IPs /U IPiso ∆fusH¯Ratio Structure (K) — — — (kJ/mol)A100 0 CsCl 1250 1.08 0.78 1.72 4406 CsCl 1200 1.08 0.81 1.57 4010 CsCl 1150 1.09 0.83 1.47 3514 CsCl 1000 1.11 0.89 1.37 2718 CsCl* 800 1.13 0.95 1.29 1922 var-CsCl 650 1.14 0.97 1.22 22A133 0 CsCl 1150 1.10 0.84 1.71 3406 CsCl 1100 1.10 0.86 1.56 3110 CsCl 1000 1.12 0.90 1.46 2614 CsCl* 850 1.13 0.95 1.37 1818 111n 700 1.08 1.00 1.30 1822 111n 650 1.10 0.98 1.22 2126 111n 550 1.10 0.99 1.18 16A167 0 NaCl 1200 1.03 0.86 1.70 3606 NaCl 1150 1.04 0.87 1.55 3310 NaCl 1100 1.05 0.89 1.44 2914 NaCl 950 1.07 0.94 1.36 2118 NaCl* 750 1.10 0.99 1.28 1822 111n 700 1.09 0.96 1.20 1826 111n 550 1.11 0.99 1.17 15Continued on next page. . .90Table 4.1: (Continued)Size d Crystal Ts ρ¯∗(s) ρ¯∗(l) U¯ IPs /U IPiso ∆fusH¯Ratio Structure (K) — — - (kJ/mol)A200 0 NaCl 1250 1.10 0.86 1.69 4006 NaCl 1200 1.11 0.88 1.54 3710 NaCl 1150 1.11 0.89 1.43 3214 NaCl 1000 1.12 0.93 1.34 2418 NaCl 800 1.13 0.98 1.26 1822 NaCl 700 1.11 0.97 1.19 1626 NaCl 600 1.11 0.96 1.16 14A233 0 NaCl 1250 1.14 0.88 1.66 3906 NaCl 1200 1.15 0.89 1.52 3710 NaCl 1100 1.16 0.92 1.42 3114 NaCl 950 1.17 0.96 1.33 2318 NaCl 800 1.14 0.98 1.25 1622 NaCl 650 1.15 0.99 1.19 1626 NaCl 550 1.13 0.99 1.15 13A267 0 NaCl 1150 1.19 0.92 1.65 3506 NaCl 1100 1.19 0.93 1.50 3210 NaCl 1050 1.19 0.94 1.40 2714 NaCl 900 1.18 0.97 1.32 1918 NaCl 750 1.15 0.99 1.24 1322 fcc(+) 600 1.16 1.00 1.18 1226 fcc(PC) 550 1.09 0.98 1.14 8A300 0 NaCl 1100 1.20 0.94 1.62 3006 NaCl 1050 1.21 0.95 1.48 2710 NaCl 900 1.22 1.00 1.39 2114 NaCl/fcc(+) 750 1.22 1.03 1.31 1518 fcc(+) 650 1.17 1.03 1.23 1022 fcc(+) 550 1.15 1.02 1.17 826 fcc(PC) 550 1.07 0.98 1.13 691The melting point trends across Set A are summarized in Fig. 4.4, wherewe show Ts to Tl for each salt. The d = 0 series of salts are the 1C - 1C salts,repeated from Fig. 4.2. The distance between counterion centers is held constant inSet A. The changes in the melting points of the 1C - 1C salts are due to the ion sizeratios and, more specifically, to the changes in the interactions between like ions. Asthe size ratio increases, the electrostatic interactions become less repulsive betweenlarger ions and more repulsive between smaller ions.For the Set A salts with d ≤ 18, the melting point trends are similar. TheA100 and A200 salts are the two optimal size ratios for the CsCl and NaCl structures,and also have the highest melting points for d ≤ 18. Each series of salts with equald ≤ 18 shows a decrease in melting point from A100 to A133 (destabilizing CsClsolid), an increase from A133 to A200 (stabilizing NaCl solid), and a decrease at orbeyond A233 (destabilizing NaCl solid). Each size ratio series shows two distinctmelting point peaks corresponding to the CsCl and NaCl solid structures whend ≤ 18.The melting point trends change for the d = 22 series of salts. There is onlyone peak in the melting temperature between A167 and A200. The A100 2L67-22- 1C salt is not stable as a CsCl solid like the A100 salts with d ≤ 18. Instead, thecrystal structure of the A100 2L67-22 - 1C salt is orientationally ordered.The variation of melting temperatures across size ratios ranges from 150 Kfor the 1C - 1C salts, to 250 K at d = 10 and 14, and down to 50 K for d = 26.For all of the size ratios except A200, Ts = 550 K. Ts of the A200 2L67-26 - 1C saltis only 50 K higher. The consistency of Ts for the d = 26 salts suggests that thevariation in melting temperature due to the size ratios (with constant σ±) can besuppressed by large cation displacement distances (d).92100 133 167 200 233 267 300Cation Size (% of Anion Size)400600800100012001400Temperature (K)d = 0d = 06d = 10d = 14d = 18d = 22d = 26Set A:   T   to T CsCl NaClorientationallyordered crystalstransitionspremeltingtransitionss-ss lFigure 4.4: Melting temperature summary for the salts in Set A. The blocks give the50 K range of Ts to Tl for each salt. The dashed grey lines separate the salts by theirsolid types. The majority freeze as NaCl solids. On the left, the salts with equal ornear-equal size ions freeze as CsCl solids, show s-s transitions, or in the bottom left,become trapped in glassy states upon cooling, depending on d. In the bottom right,the horizontal lines indicate the approximate onset of fast ion conductor (dashedlines) and plastic crystal (dotted lines) phases. (See text for description and onsetcriteria details.)93A plot of Ts against the cation displacement distance, d, for the salts in SetA is shown in Fig. 4.5, and show similar trends for all of the size ratios. The plotin Fig. 4.5, highlights that the A133 and A300 salts have the lowest melting pointswhile the A100, A200 and A233 are consistently the size ratios with the highestmelting temperatures.0 4 8 12 16 20 24d (0.01 nm)60080010001200TsA100A133A167A200A233A267A300Size Set AFigure 4.5: Melting points of the Set A salts as a function of the cation charge dis-placement distance d. The A100, A200, and A233 salts show the smallest meltingpoint reductions, while the A133 and A300 show the largest melting point reduc-tions.The melting temperature trends of the Set A salts, with constrained σ±, havebeen established and discussed. We now move on to Set B, where σ− is fixed.944.3 Size Set B (Unconstrained σ±)In Set B, we examine the impact of ion size on the melting point where σ−is held constant at 0.50 nm while σ+ and σ± increase. Because σ± increases, theattractive electrostatic interactions between counterions become weaker. As in SetA, we consider size ratios between 1:1 and 3:1, with the specific size ratios given inTable 2.1.In contrast with the charge-centered (1C - 1C) salts from Set A, the Set B1C - 1C salts do not exhibit obvious markers of the crystal structure changes intheir melting points trends. Instead, the melting points of the Set B 1C - 1C saltsdecrease monotonically. As the size ratio increases from 1:1 to 3:1, Ts of the 1C- 1C salts in Set B drops from 1250 K to 750 K, which is a significantly largerreduction than the salts in Set A (150 K). The melting point reduction of the Set B1C - 1C salts is attributed to the increasing size of the cation and resulting weakerelectrostatic interactions.Properties of the Set B salts are given in Table 4.2. The properties in Ta-ble 4.2 are the same as those listed for Set A in Table 4.1, but with an additionalcolumn for the average potential energy of the solid U¯s at Ts. The average potentialenergies of the solids U¯s at Ts, for the Set B 1C - 1C salts are listed in Table 4.2 asthe d = 0 entry, and increase from -494 to -250 kJ/mol of ion pairs as the size ratioincreases. The increase of about 245 kJ/mol of ion pairs contrasts sharply with the1C - 1C salts in Set A, which show an increase of about 30 kJ/mol of ion pairs, from-494 to -465 kJ/mol of ion pairs. The average potential energies of the Set B saltsare in the same range as some prototypic ionic liquids (between -300 to -400 kJ/molof ion pairs).11,16,7195The trends across the Set B 2L67 - 1C salts properties given in Table 4.2 arealso worth noting. The average reduced densities of the solids vary between 1.0 and1.1, much like the salts in Set A. The trend of increasing solid densities reverses forthe B267 and B300 salts with large d. These salts also show premelting transitions,like the salts in Set A discussed in Chapter 4.4. The decrease in the average re-duced density of the solid and ∆fusH¯ ≤ 10 kJ/mol of ion pairs are consistent withpremelting transitions to higher entropy states. Within each size ratio, the trendsfor U¯ IPs /U IPiso and ρ¯∗(l) are consistent with the Set A salts; as d increases, U¯ IPs /U IPisogenerally decreases while ρ¯∗(l) generally increases.Table 4.2: Properties of the Set B salts. Each salt is identified by the size ratioand cation charge displacement distance d. As in Table 4.1, the crystal structureis that of the last stable solid before melting, where fcc(+) means a fcc latticeof cations. The highest temperature where the solid persists in the two-phasesimulations is Ts. The average reduced densities, ρ¯∗, average potential energy ofthe solid, U¯s, and ∆fusH¯ are given at Ts. At each size ratio, the d = 0 entry refersto the charge-centered 1C - 1C salt.Size d Crystal Ts ρ¯∗(s) ρ¯∗(l) U¯s U¯ IPs /U IPiso ∆fusH¯Ratio Structure (K) — — kJ/mol — (kJ/mol)B133 0 CsCl 1050 1.08 0.83 -425 1.74 326 CsCl 1050 1.07 0.83 -426 1.61 3110 CsCl 1000 1.08 0.85 -430 1.52 2814 CsCl 900 1.09 0.89 -439 1.44 2318 CsCl 750 1.11 0.95 -454 1.36 1522 var-111n 600 1.12 1.00 -486 1.31 1826 var-111n 600 1.08 0.97 -524 1.24 19B167 0 NaCl 1000 0.99 0.83 -371 1.75 296 NaCl 950 1.01 0.86 -374 1.64 2710 NaCl 950 1.00 0.86 -375 1.56 2614 NaCl 900 1.01 0.88 -380 1.48 2218 NaCl 800 1.03 0.92 -389 1.42 18Continued on next page. . .96Table 4.2: (Continued)Size d Crystal Ts ρ¯∗(s) ρ¯∗(l) U¯s U¯ IPs /U IPiso ∆fusH¯Ratio Structure (K) — — kJ/mol - (kJ/mol)22 NaCl 700 1.05 0.95 -402 1.35 1426 111n 600 1.08 0.98 -424 1.30 1430 111n 550 1.09 0.98 -457 1.25 15B200 0 NaCl 950 1.05 0.84 -334 1.78 316 NaCl 950 1.05 0.84 -334 1.67 3010 NaCl 900 1.06 0.87 -337 1.61 2814 NaCl 900 1.06 0.86 -338 1.54 2618 NaCl 850 1.07 0.89 -343 1.47 2322 NaCl 800 1.07 0.90 -348 1.40 2026 NaCl 700 1.08 0.94 -359 1.34 1530 NaCl 600 1.08 0.97 -377 1.29 1334 NaCl 550 1.08 0.97 -402 1.24 12B233 0 NaCl 900 1.08 0.85 -301 1.78 296 NaCl 900 1.08 0.85 -301 1.69 2910 NaCl 900 1.08 0.85 -301 1.63 2814 NaCl 850 1.09 0.87 -304 1.57 2618 NaCl 800 1.10 0.89 -308 1.52 2422 NaCl 800 1.09 0.89 -310 1.45 2126 NaCl 700 1.11 0.93 -318 1.39 1730 NaCl 650 1.10 0.94 -325 1.33 1434 NaCl 600 1.08 0.95 -338 1.28 12B267 0 NaCl 800 1.11 0.88 -274 1.80 256 NaCl 800 1.11 0.88 -275 1.72 2510 NaCl 800 1.11 0.88 -275 1.66 2514 NaCl 800 1.11 0.88 -275 1.59 2318 NaCl 750 1.12 0.90 -278 1.55 2222 NaCl 750 1.11 0.90 -280 1.48 2026 NaCl 700 1.11 0.92 -284 1.43 1730 NaCl 650 1.11 0.94 -289 1.38 14Continued on next page. . .97Table 4.2: (Continued)Size d Crystal Ts ρ¯∗(s) ρ¯∗(l) U¯s U¯ IPs /U IPiso ∆fusH¯Ratio Structure (K) — — kJ/mol - (kJ/mol)34 NaCl 600 1.10 0.95 -297 1.32 1142 fcc(+) 500 1.08 0.97 -325 1.36 850 fcc(+) 450 1.04 0.97 -382 1.46 5B300 0 NaCl 750 1.11 0.89 -250 1.79 226 NaCl 750 1.11 0.89 -250 1.71 2110 NaCl 750 1.11 0.89 -250 1.66 2114 NaCl 700 1.13 0.91 -253 1.62 2018 NaCl 700 1.12 0.91 -254 1.57 1922 NaCl 700 1.11 0.91 -254 1.50 1726 NaCl 650 1.12 0.93 -258 1.46 1530 NaCl 650 1.09 0.93 -259 1.40 1234 NaCl 600 1.09 0.95 -265 1.36 1042 fcc(+) 550 1.02 0.95 -278 1.26 550 fcc(+) 450 1.04 0.97 -315 1.21 458 fcc(+) 450 1.02 0.95 -381 1.14 4We plot Ts to Tl for the Set B salts in Fig. 4.6. The melting point trends inSet B differ from those in Set A. The differences between Set A (Fig. 4.4) and SetB (Fig. 4.6) highlight how isolating the effect of size asymmetry on melting point isdifficult because of the different length scales involved.The melting points of the Set B charge-centered salts decrease monotonicallyas the size ratio increases. As expected, moving 2/3 of the unit charge away fromthe cation center generally reduces the melting point. At small d values, the meltingpoint trends follow the 1C - 1C trends closely across increasing size ratios.98100 133 167 200 233 267 300Cation Size (% of Anion Size)20040060080010001200Temperature (K)d = 0d = 06d = 10d = 14d = 18d = 22d = 26d = 30d = 34d = 42d = 50d = 58Set B:  T  to TCsClNaCls-stransitionsorientationallyordered crystalspremeltingtransitionss lFigure 4.6: Melting point summary of the Set B salts. The blocks give the 50 Krange of Ts to Tl for each salt. The dashed grey lines separate the salts by their solidtypes and the boundaries are similar to those found in Set A (Fig. 4.4). In the bottomleft, the salts become trapped in glassy states upon cooling. In the bottom right,the horizontal lines indicate the approximate onset of fast ion conductor (dashedlines) and plastic crystal (dotted lines) phases.99At d = 14, we see the first indication that the crystal structures affects themelting point, with the B133 2L67-14 - 1C (CsCl) and B167 2L67-14 - 1C (NaCl)salts have the same melting point ranges. As d increases to 18, the melting point ofthe B133 salt (CsCl) is lower than either of the corresponding A100 (CsCl) or B167(NaCl) salts. The drop in melting point of the B133 2L67-18 - 1C salt changes themelting point trends across the d = 18 series. The series of d = 18 melting pointsdevelops a second peak at B200 that is higher than the CsCl peak for the A1002L67-18 - 1C salt. At d ≥ 22, the second melting point peak extends to even highersize ratios, with Ts = 800 K for the size ratios of 2:1 and 2.33:1. In the d = 26,30, and 34 salt series, the melting points have a single peak that spans multiple sizeratios. The peaks extend to larger size ratios as d increases.Insight into why the melting point peak shifts to larger size ratios is providedby the ratio of the average potential energy per ion pair in the solid to the minimumpotential energy of an isolated ion pair, U¯ IPs /U IPiso. In Set A, U¯ IPs /U IPiso decreases from1.72 to 1.62 as the size ratio increased for the 1C - 1C salts. However, in Set B,U¯ IPs /U IPiso peaks at B267 with a value of 1.80 for the 1C - 1C salts. Based on thisratio, the Set B salts with larger size ratios have more thermodynamic incentive tocrystallize.If we follow the melting point trends for d = 22, 26, or 30 in Fig. 4.6, we cansee the transition between two main influences on melting point. In the lower leftof Fig. 4.6, the cation charge distribution has the largest influence on the meltingpoints. The charge distribution leads to strong ion pairing. The crystal structuresare orientationally ordered to maintain the ion pairs. If the salt is dominated byion pairing, then the melting point is low. As the cation grows larger in Set B, theoff-center charge site recedes from the cation surface and the cation charge distri-100bution becomes a weaker influence on melting point; the melting points increase.The cations in these salts are able to reorient about the lattice sites and interactfavourably with more than one nearest neighbour, which stabilizes the solid. Oncethe cation is large enough that d is a minor perturbation, the melting points of thesalts join the downward trend due to the ion sizes.The cations increase in size as the size ratio increases in the Set B salts.In the size-symmetric salts (A100), d = 22 (0.22 nm) leads to strong ion pairingbecause the cation radius is near 0.25 nm, but d = 22 does not lead to ion pairingfor B300, where the cation radius is near 0.75 nm. A plot that takes these differencesinto consideration is shown in Fig. 4.7, where Ts/Ts(1C - 1C) is plotted against thenormalized charge arm of the cation, as defined in Eq. (2.2). The trends in Fig. 4.7are similar across size ratios, and also remarkably similar to the trends in Chapter 3(Fig. 3.9) where different fractions of the cation charge are moved off-center in size-symmetric salts.1010 0.1 0.2 0.3 0.4 0.5 0.6Normalized Charge Arm Lc0. / Ts(1C - 1C)A100B133B167B200B233B267B300Size Set BFigure 4.7: Melting points (Ts) scaled by Ts of the 1C - 1C salt (of the same sizeratio) plotted against the normalized charge arm, Lc. The trends shown here aresimilar to the trends established in Chapter 3, where the fraction of off-center chargeis varied. (See Fig. 3.9.)1024.4 Solid Structure and Ion DynamicsSix of the A267 and A300 salts show premelting transitions noted in Fig. 4.4.The A267 2L67 - 1C salts with d ≥ 22 and the A300 2L67 - 1C salts with d ≥ 14display intermediate phases between the low-enthalpy solids and the liquids. In thissection, we explore the ion dynamics for the A300 series of salts, with emphasis onthe A300 2L67-22 - 1C and A300 2L67-26 - 1C salts. We compare the salts at T =550 K, which is Ts for both salts.The mean squared displacement (msd) provides insight into the particle dy-namics and is calculated as119〈∆r(t)2〉 = 〈|r(t)− r(0)|2〉, (4.1)where r(t) is the particle position at time t, and the angular brackets denote anensemble average. A log-log plot of the msd of both ions in the A300 series of saltsat 550 K is shown in Fig. 4.8. The Einstein relation119 is used to calculate thediffusion coefficients D from the msd as D = limt→∞〈|r(t)− r(0)|2〉/6t.For each of the A300 salts, the smaller anions have a higher msd than thelarger cations, as expected. All of the salts are solids at T = 550 K, and the majorityof the ions show no net translational motion. The deviations begin with d = 14,where there is a small increase in the anion msd over 0.5 - 1.0 ns time scales. For theA300 salts with d = 18 and 22, the anions show linear (diffusive) behaviour whereasthe cations do not. For these two salts, the anions move through a cation latticewith diffusion coefficients of D− = 1.548 (± 0.108) ×10−10 m2/s at T = 550 K ford = 18, and D− = 6.719 (± 0.286) ×10−10 m2/s for d = 22.1031 10 100 1000Time (ps)0.010.1110MSD (nm2)d = 0d = 06d = 10d = 14d = 18d = 22d = 26A300T = 550 KHeating NaCl  or  111n crystalsFigure 4.8: The mean square displacements for the Set A salts with a 3:1 size ratioat T = 550 K. Cation msd are marked with dashed lines and anions with solid lines.The A300 salts with d ≤ 10 are crystalline, and the msd plots are from heating NaClcrystals. For the d ≥ 14 salts, we show the msd from heating the 111n structure.For the salts with d = 18 and 22, the anions show diffusive motion while the cationsdo not. The msd of the ions in the d = 26 salt shows no net translational motion,but the anion msd is an order of magnitude larger than that of the cation. This saltexists in a plastic crystal phase at T = 550 K.104We have not attempted to determine the premelting transition temperatureaccurately here because we are interested in the s-l transition. We approximatethe onset of the fast ion conductor phase as the lowest temperature where we seean increase in the average potential energy and an anion diffusion coefficient abovean arbitrary threshold of 1 × 10−10 m2/s during the heating of the crystals. Anestimation using the Nernst-Einstein relation puts the ion conductivities in the rangeof 0.01 S/cm, a generally accepted threshold for superionic conductor behaviour.122For comparison, the α-phase of AgI, which is stable between T = 420 K and831 K, is a common example of a superionic conductor. The iodide ions form a bcclattice and the silver ions move through it. From molecular dynamics simulationsat 560 K, the diffusion coefficient of the more mobile Ag ions has been reported as3.5 ×10−9 m2/s.123 The ion diameters used in the study were 0.063 nm for silverand 0.22 nm for iodide, giving a large:small ion size ratio of 3.49.123In the A300 2L67-26 - 1C salt, which has the largest size ratio and chargedisplacement considered here, the anion and cation both have diffusion coefficientsof zero, but the msd of the anions is an order of magnitude larger than that of thecations (0.351 nm2 vs 0.0223 nm2) at T = 550 K. Here, the salt forms a plasticcrystal phase where the ions rotate as pairs; the larger cations vibrate and reorientabout lattice positions. As the cation reorients, the anion traces a path on the cationsurface and follows the movements of the off-center charge site. The anion tracesout a larger path than the cation, which explains why the msd of the anion is anorder of magnitude larger.The approximate onset temperature of the plastic crystal phase, as with thefast ion conductor phase, is determined by an increase in the average potentialenergy during the heating of the crystal, and where the msd of the anion reaches105(a) A300 1C - 1C salt(b) A300 2L67-22 - 1C salt (c) A300 2L67-26 - 1C saltFigure 4.9: Anion motion in three salts with a size ratio of 3:1. The grey spheres(cation centers) and simulation box boundaries are shown for one instant in thetrajectory. The sets of coloured points mark the positions of eight different anionson 1.0 ps intervals over the 4.0 ns trajectory. Each of these configurations are at550 K, which is Ts for the salts shown in (b) and (c).106a constant value within about 100 ps and is significantly larger than that of thecation.The different ion dynamics between the crystalline phase, fast ion conductorphase, and plastic crystal phase for three A300 salts are shown in Fig. 4.9. Thepictures show traces of the coordinates of eight anions over the 4.0 ns simulationsat T = 550 K. For all three salts, the fcc lattice of cation centers (grey spheres)and the simulation cell boundaries are drawn at one instant in the trajectory. Thepositions of eight anions are shown as points in different colours over the entire 4.0ns trajectory, and capture the anion paths through the solid. We have not shownthe off-center cation sites in (b) and (c) for simplicity. The figures were generatedfrom heating prepared crystals (a) NaCl, (b) 111n, and (c) 111n, where the cationsin (b) and (c) have rearranged into an fcc lattice.At 550 K, the A300 1C - 1C salt is well below its melting temperature (Ts= 1100 K) and is a NaCl solid. Over the course of the simulation trajectory, theanions vibrate around their lattice positions, as shown in Fig. 4.9(a). The simulationcell is shown on an angle otherwise the anion positions would be obscured by thecations. The salts shown in (b) and (c) show very different types of anion motion andpremelting behaviour. The cations only differ by the charge displacement distance,and Ts for both salts is 550 K. In Fig. 4.9(b), the A300 2L67-22 - 1C salt, theanions are moving through the cation lattice and are not localized. In contrast toFig. 4.9(b), the A300 2L67-26 - 1C salt in (c) shows that the anions do not travelthrough the cation lattice; instead, they remain localized around a single cation.The intermediate phase between the solid and liquid for this salt is characterized byrotating ion pairs (plastic crystals).107The A267 and A300 series of salts are insulators when d is small (crystalline)and very large (plastic crystals). The salts with large size ratios and d that are notbound tightly enough to become plastic crystals become fast ion conductors.The radial distribution functions (rdfs) for the A300 series of salts at T =550 K are shown in Fig. 4.10. The fcc cation lattice is maintained for all of the salts,as shown in the CM - CM rdfs in Fig. 4.10(b). The peaks become broader and shiftto the right as d increases because the salts are closer to their melting points. Thesets of rdfs involving the anion center, AM, transform from solid-like to liquid-likeas d increases. The first peaks in the CM - AM and CC - AM rdf series shift toshorter distances as d increases, which is consistent with development of directionalion pairs, as noted earlier in Chapter 3. For the fast ion conductors (d = 18 and 22),the first peaks at short distances in the CC - AM rdfs [Fig. 4.10(d)] indicate that onaverage, the anions are paired with cations. However, the time-averaged rdfs do notprovide any insights into the longevity/lifetime of any particular cation-anion pair.Likewise, at d = 26, the rdfs offer no insight that the same cation-anion pairs areleft intact for the entire simulation and are rotating as ion pairs. The ion dynamicsin the fast ion conductor and plastic crystal phases are an area worthy of furtherstudy.1080.5 1 1.5 2r (nm)0246g(r)0246g(r)d = 0d = 06d = 10d = 14d = 18d = 22d = 260246g(r)0246g(r)T = 550 K(a)  CM - AM(c) AM - AM(d) CC - AMA300(b) CM - CMFigure 4.10: Radial distribution functions for the A300 series of salts. In (a-c),the rdfs for small d are for a NaCl structure. As d increases, the rdfs involving theanion center (AM) flatten while the cation centers remain on the fcc lattice. In (d),the first CC - AM peak shifts to shorter distances and becomes more intense as dincreases. For d = 22 and 26, the time-average structural properties cannot revealthat the ion dynamics in these two salts are very different.1094.5 Chapter SummaryWe have established the melting point trends of a large set of coarse grainmodel salts considering the effects of the ion size ratio, absolute ion sizes, and cationcharge distribution. When the σ± distance is constrained (Set A) and d is small,the melting point trends reflect changes in the crystal structures. At large d, themelting points are nearly constant and show little influence from crystal structureor size ratio. When σ− is fixed and σ+ increases (Set B), the melting point trendsof the 1C - 1C salts decrease monotonically. The trends across the 2L67 - 1C saltsin Set B highlight the conversion between the two influences we consider, the ionsizes and directional ion pairing.In addition to identifying melting point trends, we find rich solid-phase be-haviours that applies for both sets of salts. Across all size ratios, cations with smalld are able to reorient with ease in the crystal. At low size ratios and small cationcharge displacements, the salts freeze as CsCl solids. At low size ratios with 2/3 ofthe cation charge far enough off-center, the salts become trapped in glassy stateswhen cooled and we find that the underlying crystal structure is orientationally or-dered. As we increase the size ratio, the salts freeze as NaCl solids. At large sizeratios and large d, the salts display two types of premelting transitions. The saltswith large d melt from fast ion conductor phases where the smaller anions movethrough an fcc cation lattice, and the salts with the largest d considered here meltfrom a plastic crystal phases composed of ion pairs rotating on an fcc lattice.Although the solid types are consistent across Set A and B, the melting pointtrends are not. The simple models are sensitive to the different length scales andproduce different melting point trends.110Chapter 5Cation Charge Geometry: SizeSymmetric 3s - 1C Salts5.1 OverviewThe salts considered in this chapter are an elaboration on the 2L - 1C saltsstudied in Chapter 3. We redistribute the cation charge over three interaction sitesinstead of two, while the anion remains a 1C ion. The three site ion geometries aredefined in Chapter 2.1. In the first set, we consider cations with the unit chargedistributed evenly over the three interaction sites (3s33 ion geometry). In the secondset, the cations have an uneven charge distribution (3s67 ion geometry) where oneoff-center charge site (CC1) carries 2/3 of the charge and the other off-center chargesite (CC2) carries the remaining 1/3 of the charge. We only consider d = 18 for thesesalts because we are interested in isolating the influence of θ, the angle formed bythe three interaction sites, CC1-CM-CC2, rather than focusing on the displacementdistance.1115.2 Melting Point ResultsThe 3s33 - 1C salts have the cation charge evenly distributed over threeinteraction sites. All of the 3s33 - 1C salts spontaneously crystallize during thehysteresis cycle. One salt, the 3s(60)33 - 1C salt, becomes trapped in a “glassy”state upon cooling, but spontaneously crystallizes during the heating part of thehysteresis cycle. The 3s33 - 1C salts all crystallize as CsCl crystals where thecations are orientationally disordered.The hysteresis plots of six 3s - 1C salts are shown in Fig. 5.1. The first columnincludes the hysteresis plots of three 3s33 - 1C salts with θ = 90◦, 135◦, and 180◦.Increasing θ does not influence the shape of the 3s33 - 1C hysteresis loops as muchas varying d did in the 2L - 1C salt series. The value of ∆fusU¯ increases slightly andthe thermal bounds on the hysteresis (T− to T+) are consistently separated by about450 K to 500 K. The change in melting temperature ranges (Ts to Tl) increases byabout 350 K as θ increases for the 3s33 - 1C salts, from Ts = 800 K for the θ = 0◦(2L67-18 - 1C) salt to Ts = 1150 K for the θ = 150◦ salt.The 3s67 - 1C salts, with an uneven charge distribution on the cation, showmore pronounced changes with increasing θ. The second column of hysteresis plotsin Fig. 5.1 are for the 3s67 - 1C salts. Of the seven 3s67 - 1C salts, three of them(with θ ≥ 135◦) spontaneously crystallize into CsCl structures upon cooling. Theother four of the 3s67 - 1C salts (with θ ≤ 120◦) become trapped in “glassy” statesupon cooling. Two of those four salts, with θ = 105◦ and 120◦, crystallize uponheating.112-580-560-540-520-500-480-460U (kJ/mol)CoolingHeatingHeat CsCltwo-phase simulations-580-560-540-520-500-480-460U (kJ/mol)300 600 900 1200Temperature (K)-580-560-540-520-500-480-460U (kJ/mol)300 600 900 1200Temperature (K)(a) 3s(90)33 (b) 3s(90)67(c) 3s(135)33 (d) 3s(135)67(e) 3s(180)33 (f) 3s(180)67Figure 5.1: The hysteresis loops of six 3s - 1C salts. All of the salts shown crystallizeas CsCl solids except the 3s(90)67-18 - 1C salt shown in (b) which vitrifies uponcooling.113Across the rows in Fig. 5.1, the θ values are constant. In the hysteresis plotswith constant θ, the redistribution of 1/3 of the charge on the cation decreasesaverage potential energies of both phases and the melting temperatures by about300 K. The magnitude of the Ts decrease is consistent with the work in Chapter 3,where the 2L33-18 - 1C has Ts = 1100 K and the 2L67-18 - 1C salt has Ts = 800 K.For all of the 3s - 1C salts, we prepared and heated a CsCl crystal until itmelted. The CsCl structure is likely the most stable solid for these size-symmetricsalts when the cations are free to reorient, a condition that is most likely satisfiednear the melting transition. At lower temperatures, the CsCl structure might notbe the most stable, especially for the 3s67 - 1C salts with θ ≤ 135◦. For the 2L - 1Csalts studied in Chapters 3 and 4, when we heated crystal structures that were notthe most stable, the plots of the average potential energies against temperature werenot smooth until the crystal rearranged or the cations gained enough thermal energyto reorient. We see similar behaviour for U¯CsCl(T ) at low temperatures for the 3s67- 1C salts with θ ≤ 135◦. In Fig. 5.1(d), the cations in the CsCl crystal undergorotational relaxations between 300 K and 550 K, after which U¯CsCl(T ) stabilizes.The fluctuations in U¯CsCl(T ) at low temperature suggests that the CsCl structuremay be metastable when the cations lack the thermal energy to reorient.The 2L100-18 - 1C salt studied in Chapter 3 would correspond to 3s(0)67-18- 1C, and melts directly from the orientationally ordered 111n crystal structure. Itis reasonable to expect, then, that some of the salts with small θ would also haveorientationally ordered crystal structure(s). However, we have not explored thecrystal structure possibilities as rigorously as we did with the 2L - 1C salts.114The properties of the 3s - 1C salts are presented in Table 5.1. We observesimilar trends in the salt properties of both sets. As θ increases, the solid densitiesare largely unaffected whereas the densities of the liquids generally decrease. Alsoas θ increases, Ts, U¯s(Ts), and ∆fusH¯ generally increase. The 3s(180)67-18 - 1C saltis an exception to the U¯s(Ts) and ∆fusH¯ trends. The linear charge arrangementpermits both off-center cation sites to interact favourably with anion neighbours,stabilizing the orientationally ordered crystal.Table 5.1: Properties of the 3s - 1C salts. The cation geometry is specified bythe charge distribution and angle formed by the interaction sites, CC1-CM-CC2.The definitions of the rest of the properties are consistent with Tables 3.1, 4.1,and 4.2.Cation Type θ Ts ρ¯∗(s) ρ¯∗(l) U¯s ∆fusH¯(degrees) (K) — — kJ/mol (kJ/mol)3s33-18 60 950 1.115 0.909 -526 2390 1050 1.102 0.872 -516 28105 1050 1.104 0.870 -514 29120 1100 1.093 0.851 -509 31135 1100 1.095 0.850 -507 31150 1150 1.086 0.831 -505 33180 1100 1.103 0.851 -508 323s67-18 60 600 1.177 1.044 -581 1290 750 1.145 0.978 -559 16105 800 1.099 0.957 -547 14120 800 1.130 0.956 -544 17135 800 1.106 0.955 -538 15150 850 1.093 0.935 -536 18180 850 1.112 0.929 -555 33115A summary of the five transition-related temperatures for both sets of 3s -1C salts is given in Fig. 5.2. We have included the 2L67-18 - 1C and 2L100-18 -1C salts as the θ = 0◦ cases. In both sets, the influence of θ on the melting pointis small when θ ≥ 105◦. When θ ≤ 90◦, the localization of charge on one side ofthe cation affects the solid-liquid transition. When θ is small, the cation can forma directional ion pair with one anion, much like the 2L - 1C salts studied earlier.Both θ = 0 cases, the 2L67-18 - 1C and 2L100-18 - 1C salts, vitrify. The3s(60)33-18 - 1C salt spontaneously crystallizes, whereas the 3s(60)67-18 - 1C saltvitrifies. The melting point range of the θ = 0 case (2L100-18 - 1C salt) is for the111n structure. The melting point of the 3s(60)67-18 - 1C salt in a CsCl structureis 200 K lower, which may be because the salt has an orientationally ordered crystalstructure with a slightly higher melting point.The 3s67-18 - 1C salts vitrify when θ ≤ 90◦. When θ ≥ 105◦, the cationsare more likely to have favourable interactions with multiple anions at once. Thecation interacting favourably with more than one anion is likely one reason why the3s67-18 - 1C salts with θ ≥ 105 spontaneously crystallize.1160 30 60 90 120 150 180theta (degrees)400600800100012001400Temperature (K)T- to T+Tg to T+Ts  to Tl 0 30 60 90 120 150 180theta (degrees)3s67-18  - 1C3s33-18  - 1CFigure 5.2: The five transition-related temperatures for the 3s - 1C salts. The saltswith uneven cation charge distribution (3s67 - 1C) have lower melting points thanthose with even cation charge distribution (3s33 - 1C). Within each set, the meltingpoints of the 3s - 1C salts with θ ≤ 90◦ vary more than the salts with θ ≥ 105◦. Onceθ ≥ 105◦, the melting points increase slightly as the angle between the off-centercharge sites increases.1175.3 Chapter SummaryWe have shown that introducing a third interaction site on the cation thatbears 1/3 of the unit charge has the greatest influence on the melting temperaturewhen the other 2/3 of the charge is off-center and the CC1-CM-CC2 angle, θ, isless than 90◦. Once θ ≥ 105◦, varying θ has a minor influence on the melting point.At large θ, the off-center cation sites can interact favourably with multiple anionswhich reduces the tendency to form discrete ion pairs.As with the 2L - 1C salts, the crystal structures of the 3s - 1C salts is stronglydependent on whether the cations are able to reorient. Most of the 3s - 1C saltscrystallize as CsCl solids. The cations of the 3s33 - 1C salts are able to reorient andthe CsCl structure is orientationally disordered. The cations of the 3s67 - 1C saltsshow a “transition” of sorts. Cations with small angles (θ ≤ 90◦) vitrify, likely withan underlying orientationally ordered crystal structure, whereas cations with largeangles (θ ≥ 105◦) crystallize as CsCl solids. The cations with a linear arrangementof interaction sites form an orientationally ordered CsCl crystal.118Chapter 6Conclusion6.1 Analysis SummaryThe work presented in this thesis is a systematic investigation into how thephysical attributes of model salts affect the melting temperature at P = 1 bar. Weconsider how the solid and liquid phases are affected by the ion charge distribu-tions, the ion size ratios, the ion sizes, and the ion symmetries by studying the saltproperties using a series of NPT molecular dynamics simulations. The simulationapproach minimizes the hysteretic effects around the phase transition. The meltingpoint ranges (Ts to Tl) of each salt are limited to 50 K, which is sufficiently accuratefor our aim of establishing melting point trends over large sets of simple model salts.In Chapter 3, we determine how distributed cation charge reduces the meltingpoints of 2L - 1C salts. Redistributing the cation charge can reduce the melting pointby over 50%, compared to the charge-centered case. Kobrak and Sandalow’s notionof the charge arm62 is a useful descriptor for correlating the melting points of the119model salts (Fig. 3.9). Moving a fraction of the cation charge off-center reduces theenthalpy of the liquid more than that of the solid, resulting in lower melting points.At large charge arms, the ions form strong, directional ion pairs and the liquidstend to vitrify. A crystal structure that accommodates ion pairs, the orientationallyordered 111n structure shown in Fig. 3.3), has a lower enthalpy than the vitreousstates.In Chapter 4, the ion size ratios and absolute ion sizes are varied for the 2L- 1C salts. The size ratios are varied between 1:1 and 3:1 in two ways. In the firstset, Set A, the distance between the cation and anion centers is fixed (σ± = 0.50nm, Set A), while the relative sizes of the ions is varied. The melting points do notvary much for the Set A 1C - 1C salts, dropping by 150 K between a 1:1 and 3:1size ratio, although changes in the melting point trends do contain signatures of thestable crystal structure changing from CsCl to NaCl. Once part of the cation chargeis moved off-center, the melting point decreases (as in Chapter 3). The melting pointtrends of the Set A 2L - 1C salts also show signs of the underlying crystal structurechanging with size ratio. At large cation charge displacements (d = 26), the meltingpoints of the salts show almost no variation between size ratios of 1.33:1 and 3:1.Only the A200 2L67-26 - 1C salt has a slightly higher (50 K) melting point range.In the second set, Set B, the size ratios are varied by keeping the anionsize fixed while increasing the size of the cation (σ− = 0.50 nm). This case ismore practical for representing real ionic liquids. In Set B, the charge-centeredsalts show a larger decrease in the melting point trends due to weaker electrostaticinteractions. Unlike the Set A 1C - 1C series, the Set B 1C - 1C salt series does notshow signs of different crystal structures affecting the melting point. Including thecation charge displacement in Set B reveals the subtle interplay between directional120ion pairs (due to cation charge displacement) and ion sizes on the melting point.The charge-centered salts in Set B provide an upper bound on the melting points,which decreases monotonically with increasing size ratio. Once part of the cationcharge is moved off-center, the melting point trends change (no longer monotonic),revealing some signs of crystal structure dependence akin to the trends establishedin Set A. The salts with the largest d considered in each size ratio up to 2:1 haveconsistent melting point ranges, between 550 K and 600 K, as in Set A. As in SetA, the salts with large size ratios and large d show premelting transitions. Of allthe salts considered in this thesis, the B300 2L67-50 - 1C and B300 2L67-58 - 1Csalts have the lowest melting point ranges, both between 450 and 500 K, with thesalts showing premelting transitions near 350 K and 300 K, respectively.The compounding effects of size ratio and charge displacement also changethe ion dynamics in the solids. The salts with size ratios of 2.67:1 and 3:1 showpremelting transitions for salts with large d. The A267 and A300 salts with thelargest d cations become plastic crystals before melting, where ion pairs rotate aboutfcc lattice positions. The salts with slightly smaller cation charge displacementsbecome fast ion conductors, where the smaller anions move through fcc lattices ofcations. Unlike typical fast ion conductors, the phase found here is not induced byion size differences alone, but by the combination of disparate ion sizes and cationcharge displacement.In Chapter 5, the cation charge is distributed over three sites instead of two,and find that θ, the angle made by the CC1-CM-CC2 interaction sites on the cation,is a minor influence on the melting point. Two different cation charge distributionsare considered. The unit charge is evenly distributed in the series of 3s33 cations,where each interaction site carries 1/3 of the unit charge. In the 3s67 cations, the121charge is distributed unevenly; the two off-center sites carry 2/3 and 1/3 of the unitcharge and the cation center is uncharged. For both sets of 3s - 1C salts, the largestmelting point reduction occurs when θ ≤ 90◦. In the 3s33 - 1C series of salts, themelting point of the 3s(60)33-18 - 1C salt is 150 K lower than the 3s(180)33-18 - 1Csalt. In the 3s67 - 1C series of salts, the melting point of the 3s(60)67-18 - 1C saltis 250 K lower than the 3s(180)67-18 - 1C salt. In both sets of 3s - 1C salts, onceθ ≥ 105◦, the melting points increase only slightly as θ approaches 180◦.This work is a broad survey of coarse grain model salts and their meltingpoint trends.6.2 Connection with Experimental WorkWhile we are interested in isolating the impact on the melting point trendsof coarse grain model salts, it is insightful to see how well our results compare withexperimentally observed melting point trends. We compare our salts to a set ofalkyl-substituted ammonium bromide salts.We take NH4Br as the reference salt, which has a normal melting point of725 K and is somewhat analogous to the A100 1C - 1C salt. The NH+4 and Br –ions have centered (or symmetrically distributed) charges, similar sizes, and arespherical (or nearly so). The series of tetraalkylammonium bromide salts are castinto the place of the 1C - 1C salts in size Set B (in Chapter 4), where the anion isfixed (Br – ) and the cation increases in size. The length of the alkyl chain functionsanalogously to the size ratio in size Set B, where the cation size increases as thelength of R2 increases. We are assuming that three R2 alkyl tails keep the cationshape nearly spherical, which is a reasonable assumption for chain lengths of four or122fewer carbons because there is little flexibility in the alkyl conformations.124 Longeralkyl chains are unable to maintain a sphere-like ion shape due to the multitude ofconformations available and the elongation of the ion shape.To vary the charge distribution of the cation, we consider the homologousseries of R1(R2)3NBr salts, where R1 and R2 are either hydrogen atoms or n-alkylchains. The length of R1 is varied and functions somewhat analogously to varyingd. When R1 and R2 are the same length, the d = 0. When R1 is shorter than R2, dincreases (the nitrogen atom is closer to the ion surface in one direction).The melting points of the alkyl-substituted ammonium bromide salts aregiven in Table 6.1 and plotted against the length of R2 in Fig. 6.1.Table 6.1: Normal melting points of substituted ammonium bromide salts. Theseries of R1(R2)3NBr salts include size effects and charge asymmetry. The bromidesalts where the cation has n-butyl chains or shorter are comparable to the sizeSet B salts from Chapter 4. The size ratios were calculated from ionic radii.125,126The melting points that are listed as sourced from Material Safety Data Sheets(MSDS) are from Sigma Aldrich.Cation Melting Point Cation:Br – SourceR1 R2 (K) Size RatioH H 725 0.93 MSDSC2H5 H 434 127n-C3H7 H 456 127n-C7H15 H 490 127CH3 CH3 573 1.77 MSDSC2H5 CH3 603 53n-C3H7 CH3 513 53n-C4H9 CH3 468 53n-C7H15 CH3 455 55CH3 C2H5 573 53Continued on next page. . .123Table 6.1: (Continued)Cation Melting Point Cation:Br – SourceR1 R2 (K) Size RatioC2H5 C2H5 558 2.12 MSDSn-C4H9 C2H5 486 55n-C5H11 C2H5 419 53n-C7H15 C2H5 393 53n-C3H7 n-C3H7 543 2.53 MSDSn-C4H9 n-C3H7 494 MSDSCH3 n-C4H9 393 53n-C4H9 n-C4H9 393 2.75 53The melting points of the tetraalkylammonium bromide salts, shown inFig. 6.1 with black borders, decrease monotonically with increasing alkyl chainlength, which is similar to the trends observed for the series of 1C - 1C salts inSet B of Chapter 4. The cation:anion size ratios span the same region exploredin Chapter 4 where ammonium bromide and tetrabutylammonium bromide haveapproximate size ratios of 0.93 and 2.75, respectively.For the cations where R1 6=R2, the ion size, shape, mass, mass distribution,and charge distribution change with respect to the tetraalkyl substituted cation.The normal melting points of ethylammonium bromide (Tm = 434 K) and n-propylammonium bromide (Tm = 456 K) are lower than that of ammonium bro-mide (Tm = 725 K) by 250-300 K. The changes in ion size, shape, and chargedistribution lower the melting point relative to ammonium bromide. Following theR1=ethyl (green blocks), propyl (yellow blocks), and butyl (red blocks) across in-124H CH3 C2H5 n-C3H7 n-C4H9R2300400500600700800Normal Melting Temperature (K)R1 = HR1 = CH3R1 = C2H5R1 = C3H7R1 = C4H9R1 = C7H15R1 = R2Figure 6.1: Plot of the normal melting points of R1(R2)3NBr substituted ammo-nium bromide salts. The melting points are plotted as a function of the trialkylsubstituents on the cation. The cation sizes increase from left to right.125creasing R2 length in Fig. 6.1 gives a generally increasing melting point trend untilthe tetraalkylammonium bromide salts are reached. The melting points of the saltswith R1=butyl cations decrease from R2=n-propyl to the tetrabutylammonium bro-mide salt, likely due to size effects.Three points are included in Fig. 6.1 for R1=n-heptyl, which are marked withgrey crosses. The trends for these salts do not follow the trends from Chapter 4,because the ion shapes are far from spherical and the ion sizes/size ratios are outsidethe region studied.The melting point trends of the coarse grain model salts studied in thisthesis can be connected to those of a series of real salts, shown here for examplealkyl substituted ammonium bromide salts. The connection demonstrates the valueof studying coarse grain salts and that the general insights are transferable to realsalts, even though the coarse grain cations are much simpler than the ammoniumcations and the coarse grain anions are much larger than bromide ions.6.3 Future WorkThis project has fostered some ideas for new research questions that areadjacent to the current project. The work in Chapter 3 raises interesting questionsabout the rest of the phase diagram of these simple model salts. For example,how do the triple points and critical points vary with cation charge displacement?Ganzenmu¨ller and Camp have identified the critical points of similar simple modelsalts.57 Additionally, it would be interesting to determine where the s-s transitionsbetween the ordered and disordered solids occur.In earlier work,92 the dynamical properties of molten model salts under NV T126conditions at T = 1200 K were studied. It would be interesting to establish trendsfor dynamic properties like diffusion, viscosity, and electrical conductivity, at a con-sistent temperature like 100 K above each salt’s melting point. A detailed profile ofthe directional ion pairing (formation/degradation rates, longevity of ion pairs) ator just above the melting point would also be insightful. The work done on thesesimple models60,61,92,128 has reinforced the importance of directional ion pairing onthe dynamic and thermodynamic properties.The work discussed in Chapter 4 generated some interesting questions. Itwould be of interest to explore the fast ion conductor phases of the ions with large sizeratios and large cation charge displacements. One could characterize the premeltingsolid phase using MD simulations (with more frequent output) to calculate thecurrent-current correlation functions and the electrical conductivities. It would alsobe interesting to consider a series of A300 salts with smaller d increments to refinethe location of the transition from a fast-ion conductor premelting phase to a plasticcrystal premelting phase across salt models. Studying these properties might givesome further insight into features of molecular ionic salts that make them attractiveas solid state electrolytes.80Also from Chapter 4, the transition from a fast-ion conductor phase to aplastic crystal phase might offer some interesting insights into the heterogeneousdynamics found in the liquid phase of ionic liquids. We have demonstrated that theanion dynamics in the solids are sensitive to a change of 0.04 nm in the cation chargedistribution. If a cation with conformational flexibility could sample molecularconfigurations with charge displacements between these two cases [A300 salts withd = 22 (fast ion conductor) and d = 26 (plastic crystal)], then the salt would have amix of paired (localized) and unpaired (delocalized) ions, and might show dynamic127heterogeneity in the solid. And if, (and this is a big if) similar behaviour is observedin the liquid, then dynamic heterogeneities would be expected, and would act as asignature of the anions responding to different cation conformations. It would beinteresting to study the liquid dynamics of neat molten salts just above the meltingpoint, and to contrast the dynamics with cases where the single cation type wasreplaced with a binary or ternary mixture of closely-related cation types to emulatea distribution of cation conformations.The MD method used to establish the melting points was sufficient for es-tablishing the melting point trends across a set of salts, but lacks precision due toa dependence on the barostat in the two-phase simulations. There are two relatedtopics that would merit further investigation. It might be useful to calculate themelting point of the 1C - 1C salts in this work using a more rigorous free energymethod. Making a comparison with an “exact” method might give insight into howto better select barostat parameters. A second, related project would aim to quan-tify, minimize, or eliminate the melting point sensitivity to the barostat parameters,which would likely involve larger system sizes.The coarse grain models employed here reproduce some interesting aspectsof ILs. There are many possible ways to make the simple salts more realistic, andalso other ways to explore theoretically interesting model salts. In order to makethe ion models more realistic, one could• explore different ion shapes,• use a polarizable model,• consider ions with internal flexibility (multiple conformations), or• create cation-anion pairs using ion topologies that resemble well-studied ILs.128The exploration of different ion shapes is a natural extension of this work.Crystal structure determination of particles with various shapes has been taken upby Dijkstra and coworkers.114,129–134 Research on the s-l transitions of ions shaped asfused spheres, prolate/oblate ellipsoids, and discs and the properties of such shapeswhen paired with a spherical counterion would be particularly insightful. It wouldalso be interesting to add (flexible) chains of LJ spheres to the ions studied here tomimic alkyl tails of various lengths. In the present work, many different solidificationbehaviours are observed, from crystallization as “regular” crystals, crystals wherethe cation reoriented within the solid, plastic crystals, and fast ion conductors, tovitrification. Liquid crystal phases are notably absent, but the ions studied here arenot aspherical enough to form liquid crystal phases. A similar approach could beused to categorize and characterize simple models for liquid crystal ionic liquids.In terms of adapting the models for theoretical purposes, one could• systematically reduce the charge magnitudes on the ions from ±1 e to about±0.6 e. Reducing the magnitude of the ion charge (charge scaling) has beenexplored as a way of reconciling dynamical property differences between sim-ulations and experiments,12,58,65,135• vary the LJ interaction strength ε in two ways, (a) vary the magnitude and(b) explore ions with different size ratios keeping ε per volume or surface areaconstant,• vary the ion topologies by adding more interaction sites.A practical project could look at the mapping between the trends found hereand homologous series of ILs, or making comparisons with experimental work.129Bibliography[1] H. Weinga¨rtner, Angew. Chem., Int. Ed. 47, 654 (2008).[2] S. Jayaraman and E. J. Maginn, J. Chem. Phys. 127, 214504 (2007).[3] N. V. Plechkova and K. R. Seddon, Chem. Soc. Rev. 37, 123 (2008).[4] R. Hayes, S. Imberti, G. G. Warr, and R. Atkin, Phys. Chem. Chem. Phys.13, 3237 (2011).[5] A. R. Katritzky, R. Jain, A. Lomaka, R. Petrukhin, M. Karelson, A. E. Visser,and R. D. Rogers, J. Chem. Inf. Comput. Sci. 42, 225 (2002).[6] E. W. Castner and J. F. Wishart, J. Chem. Phys. 132, 120901 (2010).[7] I. Krossing, J. M. Slattery, C. Daguenet, P. J. Dyson, A. Oleinikova, andH. Weinga¨rtner, J. Am. Chem. Soc. 128, 13427 (2006).[8] B. L. Bhargava, S. Balasubramanian, and M. L. Klein, Chem. Commun. 29,3339 (2008).[9] E. J. Maginn, J. Phys.: Condens. Matter 21, 373101 (2009).[10] U. Preiss, S. Bulut, and I. Krossing, J. Phys. Chem. B 114, 11133 (2010).[11] S. Tsuzuki, H. Tokuda, K. Hayamizu, and M. Watanabe, J. Phys. Chem. B109, 16474 (2005).130[12] Y. Wang, W. Jiang, T. Yan, and G. A. Voth, Acc. Chem. Res. 40, 1193(2007).[13] A. R. Katritzky, A. Lomaka, R. Petrukhin, R. Jain, M. Karelson, A. E. Visser,and R. D. Rogers, J. Chem. Inf. Comput. Sci. 42, 71 (2002).[14] J. F. Brennecke and E. J. Maginn, AIChE J. 47, 2384 (2001).[15] M. N. Kobrak, Adv. Chem. Phys. 139, 85 (2008).[16] E. A. Turner, C. C. Pye, and R. D. Singer, J. Phys. Chem. A 107, 2277(2003).[17] K. Marsh, J. Boxall, and R. Lichtenthaler, Fluid Phase Equilib. 219, 93(2004).[18] R. D. Rogers and K. R. Seddon, Science 302, 792 (2003).[19] M. J. Earle, S. P. Katdare, and K. R. Seddon, Org. Lett. 6, 707 (2004).[20] C. Chiappe and D. Pieraccini, J. Phys. Org. Chem. 18, 275 (2005).[21] M. J. Marczewski, B. Stanje, I. Hanzu, M. Wilkening, and P. Johansson, Phys.Chem. Chem. Phys. 16, 12341 (2014).[22] C. P. Fredlake, J. M. Crosthwaite, D. G. Hert, S. N. V. K. Aki, and J. F.Brennecke, J. Chem. Eng. Data 49, 954 (2004).[23] M. E. Valkenburg, R. L. Vaughn, M. Williams, and J. S. Wilkes, Thermochim-ica Acta 425, 181 (2005).[24] A. S. Arico, P. Bruce, B. Scrosati, J.-M. Tarascon, and W. van Schalkwijk,Nat. Mater. 4 (2005).131[25] J. F. Wishart, Energy Environ. Sci. 2, 956 (2009).[26] S. Aparicio, M. Atilhan, and F. Karadas, Ind. Eng. Chem. Res. 49, 9580(2010).[27] P. Brown, C. Butts, R. Dyer, J. Eastoe, I. Grillo, F. Guittard, S. Rogers, andR. Heenan, Langmuir 27, 4563 (2011).[28] P. Brown, C. P. Butts, J. Eastoe, I. Grillo, C. James, and A. Khan, J. ColloidInterface Sci. 395, 185 (2013).[29] C. Ye, W. Liu, Y. Chen, and L. Yu, Chem. Commun. 21, 2244 (2001).[30] C.-M. Jin, C. Ye, B. S. Phillips, J. S. Zabinski, X. Liu, W. Liu, and J. M.Shreeve, J. Mater. Chem. 16, 1529 (2006).[31] M. Palacio and B. Bhushan, Tribology Letters 40, 247 (2010).[32] X. Zhang, X. Zhang, H. Dong, Z. Zhao, S. Zhang, and Y. Huang, EnergyEnviron. Sci. 5, 6668 (2012).[33] T. M. Chang, L. X. Dang, R. Devanathan, and M. Dupuis, J. Phys. Chem. A114, 12764 (2010).[34] F. Karadas, M. Atilhan, and S. Aparicio, Energy Fuels 24, 5817 (2010).[35] M. Hasib-ur Rahman, M. Siaj, and F. Larachi, Chem. Eng. Process. ProcessIntensif. 49, 313 (2010).[36] C. Wang, H. Luo, D.-e. Jiang, H. Li, and S. Dai, Angew. Chem. 122, 6114(2010).[37] J. G. Huddleston and R. D. Rogers, Chem. Commun. 16, 1765 (1998).132[38] A. E. Visser, R. P. Swatloski, W. M. Reichert, R. Mayton, S. Sheff,A. Wierzbicki, J. H. Davis, Jr., and R. D. Rogers, Chem. Commun. 1, 135(2001).[39] A. E. Visser, R. P. Swatloski, S. T. Griffin, D. H. Hartman, and R. D. Rogers,Sep. Sci. Technol. 36, 785 (2001).[40] R. Sheldon, Chem. Commun. 23, 2399 (2001).[41] C. M. Gordon, Applied Catalysis A: General 222, 101 (2001), CelebrationIssue.[42] T. Welton, Coord. Chem. Rev. 248, 2459 (2004), Vignettes of HomogeneousCatalysis.[43] V. I. Paˆrvulescu and C. Hardacre, Chem. Rev. 107, 2615 (2007).[44] H. Olivier-Bourbigou, L. Magna, and D. Morvan, Applied Catalysis A: General373, 1 (2010).[45] M. Malvaldi and C. Chiappe, J. Phys.: Condens. Matter 20, 035108 (2008).[46] A.-P. Hynninen and A. Z. Panagiotopoulos, Mol. Phys. 106, 2039 (2008).[47] A. A. H. Pa´dua, M. F. Costa Gomes, and J. N. A. Canongia Lopes, Acc.Chem. Res. 40, 1087 (2007).[48] W. Xu, E. I. Cooper, and C. A. Angell, J. Phys. Chem. B 107, 6170 (2003).[49] G. Raabe and J. Ko¨hler, J. Chem. Phys. 128, (2008).[50] K. Wendler, F. Dommert, Y. Y. Zhao, R. Berger, C. Holm, and L. Delle Site,Faraday Discuss. 154, 111 (2012).133[51] A. R. Katritzky, R. Jain, A. Lomaka, R. Petrukhin, U. Maran, and M. Karel-son, Cryst. Growth Des. 1, 261 (2001).[52] S. Trohalaki, R. Pachter, G. W. Drake, and T. Hawkins, Energy Fuels 19,279 (2005).[53] D. M. Eike, J. F. Brennecke, and E. J. Maginn, Green Chem. 5, 323 (2003).[54] I. Lo´pez-Martin, E. Burello, P. N. Davey, K. R. Seddon, and G. Rothenberg,ChemPhysChem 8, 690 (2007).[55] A. Varnek, N. Kireeva, I. V. Tetko, I. I. Baskin, and V. P. Solov’ev, J. Chem.Inf. Model. 47, 1111 (2007).[56] S. Alavi and D. L. Thompson, J. Chem. Phys. 122, (2005).[57] G. C. Ganzenmu¨ller and P. J. Camp, Condens. Matter Phys. 14, 1 (2011),Available at http://arxiv.org/abs/1202.4279.[58] D. Roy and M. Maroncelli, J. Phys. Chem. B 114, 12629 (2010).[59] D. Roy, N. Patel, S. Conte, and M. Maroncelli, J. Phys. Chem. B 114, 8410(2010).[60] H. V. Spohr and G. N. Patey, J. Chem. Phys. 130, 104506 (2009).[61] H. V. Spohr and G. N. Patey, J. Chem. Phys. 132, 154504 (2010).[62] M. Kobrak and N. Sandalow, Molten Salts XIV, Electrochemical SocietyPennington, NJ, 2006.[63] S. M. Urahata and M. C. C. Ribeiro, J. Phys. Chem. Lett. 1, 1738 (2010).134[64] O. Patsahan and A. Ciach, Phys. Rev. E 86, 031504 (2012).[65] E. I. Izgorodina, Phys. Chem. Chem. Phys. 13, 4189 (2011).[66] M. G. Del Po´polo and G. A. Voth, J. Phys. Chem. B 108, 1744 (2004).[67] S. Zahn, J. Thar, and B. Kirchner, J. Chem. Phys. 132, 124506 (2010).[68] P. Wasserscheid and T. Welton, editors, Ionic Liquids in Synthesis, Wiley,2006.[69] H. Feng, J. Zhou, and Y. Qian, J. Chem. Phys. 135, 144501 (2011).[70] J. N. A. Canongia Lopes and A. A. H. Pa´dua, J. Phys. Chem. B 110, 3330(2006).[71] S. Zahn, G. Bruns, J. Thar, and B. Kirchner, Phys. Chem. Chem. Phys. 10,6921 (2008).[72] D. MacFarlane, P. Meakin, J. Sun, N. Amini, and M. Forsyth, J. Phys. Chem.B 103, 4164 (1999).[73] L. Jin, K. M. Nairn, C. M. Forsyth, A. J. Seeber, D. R. MacFarlane, P. C.Howlett, M. Forsyth, and J. M. Pringle, J. Am. Chem. Soc. 134, 9688 (2012).[74] D. Hwang, D. Kim, S. Jo, V. Armel, D. R. MacFarlane, D. Kim, and S.-Y.Jang, Sci. Rep. (2013).[75] D. R. MacFarlane, J. Huang, and M. Forsyth, Nature 402, 792 (1999).[76] K. Matsumoto, R. Hagiwara, Z. Mazej, P. Benkicˇ, and B. Zˇemva, Solid StateSciences 8, 1250 (2006).135[77] E. W. Castner, C. J. Margulis, M. Maroncelli, and J. F. Wishart, Annu. Rev.Phys. Chem. 62, 85 (2011).[78] D. R. MacFarlane and M. Forsyth, Advanced Materials 13, 957 (2001).[79] J. M. Pringle, P. C. Howlett, D. R. MacFarlane, and M. Forsyth, J. Mater.Chem. 20, 2056 (2010).[80] F. Chen, S. W. de Leeuw, and M. Forsyth, J. Phys. Chem. Lett. 4, 4085(2013).[81] H. A. Every, A. G. Bishop, D. R. MacFarlane, G. Oradd, and M. Forsyth, J.Mater. Chem. 11, 3031 (2001).[82] C. M. Gordon, J. D. Holbrey, A. R. Kennedy, and K. R. Seddon, J. Mater.Chem. 8, 2627 (1998).[83] J. D. Holbrey and K. R. Seddon, J. Chem. Soc., Dalton Trans. , 2133 (1999).[84] J. D. Holbrey, W. M. Reichert, M. Nieuwenhuyzen, S. Johnson, K. R. Seddon,and R. D. Rogers, Chem. Commun. 14, 1636 (2003).[85] A. M. Scurto and W. Leitner, Chem. Commun. , 3681 (2006).[86] Y. Shimizu, Y. Ohte, Y. Yamamura, K. Saito, and T. Atake, J. Phys. Chem.B 110, 13970 (2006), PMID: 16836349.[87] R. D. Chirico, V. Diky, J. W. Magee, M. Frenkel, and K. N. Marsh, PureAppl. Chem. 81, 791 (2009).[88] L. Zheng, S.-N. Luo, and D. L. Thompson, J. Chem. Phys. 124, 154504(2006).136[89] Y. Zhang and E. J. Maginn, J. Chem. Phys. 136, 144116 (2012).[90] S. Wang, G. Zhang, H. Liu, and H. Song, J. Chem. Phys. 138, 134101 (2013).[91] S.-N. Luo, A. Strachan, and D. C. Swift, J. Chem. Phys. 120, 11640 (2004).[92] H. V. Spohr and G. N. Patey, J. Chem. Phys. 129, 064517 (2008).[93] Q. Yan and J. J. de Pablo, Phys. Rev. Lett. 88, 095504 (2002).[94] R. A. Kuharski and P. J. Rossky, J. Am. Chem. Soc. 106, 5794 (1984).[95] S.-N. Luo, T. J. Ahrens, T. C¸ag˘ ın, A. Strachan, W. A. Goddard, and D. C.Swift, Phys. Rev. B 68, 134206 (2003).[96] L. Zheng and D. L. Thompson, J. Chem. Phys. 125, 084505 (2006).[97] D. Van Der Spoel, E. Lindahl, B. Hess, G. Groenhof, A. E. Mark, and H. J. C.Berendsen, J. Comput. Chem. 26, 1701 (2005).[98] B. Hess, C. Kutzner, D. van der Spoel, and E. Lindahl, J. Chem. TheoryComput. 4, 435 (2008).[99] S. Pronk, S. Pll, R. Schulz, P. Larsson, P. Bjelkmar, R. Apostolov, M. R.Shirts, J. C. Smith, P. M. Kasson, D. van der Spoel, B. Hess, and E. Lindahl,Bioinformatics 29, 845 (2013).[100] S. Nose´, Mol. Phys. 52, 255 (1984).[101] S. Nose´, J. Chem. Phys. 81, 511 (1984).[102] W. G. Hoover, Phys. Rev. A 31, 1695 (1985).[103] M. Parrinello and A. Rahman, J. Appl. Phys. 52, 7182 (1981).137[104] S. Nose´ and M. Klein, Mol. Phys. 50, 1055 (1983).[105] U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee, and L. G. Ped-ersen, J. Chem. Phys. 103, 8577 (1995).[106] B. Hess, H. Bekker, H. J. C. Berendsen, and J. G. E. M. Fraaije, J. Comput.Chem. 18, 1463 (1997).[107] E. G. Noya, C. Vega, and E. de Miguel, J. Chem. Phys. 128, 154507 (2008).[108] C. Hardacre, J. D. Holbrey, M. Nieuwenhuyzen, and T. G. A. Youngs, Acc.Chem. Res. 40, 1146 (2007).[109] C. Vega, J. L. F. Abascal, C. McBride, and F. Bresme, J. Chem. Phys. 119,964 (2003).[110] G. C. McNeil-Watson and N. B. Wilding, J. Chem. Phys. 124, 064504 (2006).[111] A. Criado, M. Jime´nez-Ruiz, C. Cabrillo, F. J. Bermejo, R. Ferna´ndez-Perea,H. E. Fischer, and F. R. Trouw, Phys. Rev. B 61, 12082 (2000).[112] R. Rey, J. Phys. Chem. B 112, 344 (2008).[113] H. T. Stokes and D. M. Hatch, J. Appl. Crystallogr. 38, 237 (2005).[114] L. Filion and M. Dijkstra, Phys. Rev. E 79, 046714 (2009).[115] S. Yamanaka, T. Yasunaga, K. Yamaguchi, and M. Tagawa, J. Mater. Chem.19, 2573 (2009).[116] S. Yamanaka, K. Umemoto, Z. Zheng, Y. Suzuki, H. Matsui, N. Toyota, andK. Inumaru, J. Mater. Chem. 22, 10752 (2012).138[117] M. Bykov, E. Bykova, S. van Smaalen, L. Dubrovinsky, C. McCammon,V. Prakapenka, and H.-P. Liermann, Phys. Rev. B 88, 014110 (2013).[118] J. Hansen and I. McDonald, Theory of Simple Liquids: with Applications toSoft Matter, Elsevier Science, 2013.[119] M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids, Clarendonpress Oxford, 1987.[120] D. Frenkel and B. Smit, Understanding Molecular Simulation: From Algo-rithms to Applications, Computational Science Series, Elsevier Science, 2001.[121] M. Zuriaga, M. Carignano, and P. Serra, J. Chem. Phys. 135, 044504 (2011).[122] D. A. Keen, J. Phys.: Condens. Matter 14, R819 (2002).[123] A. Fukumoto, A. Ueda, and Y. Hiwatari, J. Phys. Soc. Jpn. 51, 3966 (1982).[124] H. Jin, B. O’Hare, J. Dong, S. Arzhantsev, G. A. Baker, J. F. Wishart, A. J.Benesi, and M. Maroncelli, J. Phys. Chem. B 112, 81 (2008), PMID: 18069817.[125] D. H. Aue, H. M. Webb, and M. T. Bowers, J. Am. Chem. Soc. 98, 318(1976).[126] B. Trend, D. Knoll, M. Ueno, D. Evans, and V. Bloomfield, Biophys. J. 57,829 (1990).[127] J. Tsau and D. Gilson, J. Phys. Chem. 72, 4082 (1968).[128] E. K. Lindenberg and G. N. Patey, J. Chem. Phys. 140, (2014).[129] A.-P. Hynninen and M. Dijkstra, Phys. Rev. Lett. 94, 138303 (2005).139[130] M. E. Leunissen, C. G. Christova, A.-P. Hynninen, C. P. Royall, A. I. Camp-bell, A. Imhof, M. Dijkstra, R. van Roij, and A. van Blaaderen, Nature 437,235 (2005).[131] A. P. Gantapara, J. de Graaf, R. van Roij, and M. Dijkstra, Phys. Rev. Lett.111, 015501 (2013).[132] K. Milinkovic´, M. Dennison, and M. Dijkstra, Phys. Rev. E 87, 032128 (2013).[133] M. Marechal, A. Patti, M. Dennison, and M. Dijkstra, Phys. Rev. Lett. 108,206101 (2012).[134] M. Dennison, K. Milinkovic´, and M. Dijkstra, J. Chem. Phys. 137, (2012).[135] M. Kla¨hn, A. Seduraman, and P. Wu, J. Phys. Chem. B 112, 10989 (2008).[136] D. M. Heyes, Phys. Rev. B 49, 755 (1994).[137] N. Galamba, C. A. Nieto de Castro, and J. F. Ely, J. Phys. Chem. B 108,3658 (2004).[138] H. Lee and W. Cai, Ewald summation for Coulomb interactions in a periodicsupercell, 2009, micro.stanford.edu/mediawiki/images/4/46/Ewald notes.pdf.[139] I. Souza and J. Martins, Phys. Rev. B 55, 8733 (1997).[140] S. Forrest, Science 261, 872 (1993).[141] D. M. Deaven and K. M. Ho, Phys. Rev. Lett. 75, 288 (1995).[142] D. Deaven, N. Tit, J. Morris, and K. Ho, Chem. Phys. Lett. 256, 195 (1996).[143] C. W. Glass, A. R. Oganov, and N. Hansen, Comput. Phys. Commun. 175,713 (2006).140[144] J. Niesse and H. R. Mayne, Chem. Phys. Lett. 261, 576 (1996).[145] J. A. Niesse and H. R. Mayne, J. Chem. Phys. 105, 4700 (1996).[146] A. R. Oganov and C. W. Glass, J. Chem. Phys. 124, (2006).[147] A. R. Oganov and C. W. Glass, J. Phys.: Condens. Matter 20, 064210 (2008).[148] D. C. Lonie and E. Zurek, Comput. Phys. Commun. 182, 372 (2011).[149] T. A. Weber and F. H. Stillinger, Phys. Rev. B 32, 5402 (1985).[150] W. Kob and H. C. Andersen, Phys. Rev. E 51, 4626 (1995).[151] J. R. Ferna´ndez and P. Harrowell, Phys. Rev. E 67, 011403 (2003).141Appendix AConfigurational CalculationsThere are elements that both programs written for this project, a moleculardynamics (MD) program and a genetic algorithm (GA), have in common. Generaloverviews of the program structures are shown in Fig. A.1. The calculations thatapply for both the GA and MD programs are described in this appendix. Specifi-cally, the parameters and representations of the container, the details of the particlerepresentations, and finally the techniques used to calculate the energies, forces andpressures for a single configuration are presented. These elements are consistentacross both programs written for this project and are discussed together here, andare shown in purple boxes near the center of Fig. A.1.For both programs written for this project, the energy evaluation routinesused Cartesian coordinates for the particles, whereas the propagation routines (throughconfiguration space and time) were carried out using fractional coordinates.Details specific to the MD program are presented in Appendix B. Some ofthe methods used in the MD code produced for this project as well as the techniques142implemented with the Gromacs software package are noted.Appendix C presents the specifics of the genetic algorithm and includes de-tails about the initialization, propagation, and termination of the simulation.A.1 Simulation CellThe simulation cell is a parallelepiped, with the size and shape defined bythree lengths and three angles. The lengths of the cell sides are a, b, c, and thethree angles are α, β, and γ, which are the angles between b and c, a and c, and aand b, respectively. In a Cartesian setting, the cell boundaries are usually definedby three vectors, a, b, and c, that are represented in a 3 × 3 matrix, h. The nineelements in the matrix can always be reduced to 6 non-zero elements; the matrixincludes redundant information. The conventional Cartesian setting of the cell isthat a coincides with the positive x axis, b lies in the xy plane, and c has somecomponent in the positive z direction,h =ax ay azbx by bzcx cy cz=ax 0 0bx by 0cx cy cz=|a| 0 0|b| cos γ |b| sin γ 0|c| cos β |c|(cosα−cosβ cos γ)sin γ V|a| |b| sin γ(A.1)where V is the cell volume.Another approach to reduce the 9 matrix elements to 6 and eliminate theorientational dependence of the vectors is to use the metric tensor, whose elementsare defined asg ≡ hTh. (A.2)143Molecular Dynamics Initializationinput, run parame-ters/constraints,initial configuration,N ∼ 102 − 103compute Cartesiancoordinates forconfiguration(s)for each configura-tion, calculate ener-gies, forces, pressurescorrect configuration,apply constraintsend sim-ulation?predict configu-ration at t + ∆tfinal outputstop MDyesnoGenetic Algorithm Initializationinput,run parame-ters/constraints,N ∼ 101generate a populationof unit cells (frac-tional coordinates)rank configurationsendsearch?apply evolu-tion operatorscheck configurationsfinal outputstop GAyesnoFigure A.1: MD (red) and GA (blue) Program Structures. The purple boxes arethe aspects that are common to both programs.144The metric tensor is symmetric, with the diagonal elements containing informationabout the cell lengths and the off-diagonal elements relaying information about theangles formed between those vectors. To reduce the metric tensor to the conventionalCartesian representation mentioned above, the transformation ish =√g11 0 0g11h11√g22 − h221 0g13h11(g23−h21 h31)h22√g33 − h231 − h232. (A.3)The volume of the cell can be calculated using h or g, asV = deth =√det g. (A.4)A.2 Particle RepresentationsThe position of the center of mass of each particle in a configuration is storedas a fractional value of the simulation cell axes. Fractional coordinates are denotedwith a superscript F , any other coordinates are Cartesian. A column vector ofparticle i’s center of mass fractional coordinates, rFi , is easily converted into theglobal Cartesian coordinates, ri, using the cell matrix, hT asri = hTrFi . (A.5)145The center of mass representation is all that is required for the 1C ions. For the ionswith more than one interaction site, the location of the off-center interaction sitesmust be specified. Each particle is treated as a rigid collection of interaction sites.Each interaction site within each particle has an amount of charge and particularLJ attributes assigned to it, as established in Chapter 2. The particles are rigidentities, meaning the lengths and angles between interaction sites are fixed and theinteractions between sites on the same particle are ignored. For each particle type,the location of each interaction site is defined in a particle coordinate system (SeeFig. 2.3 as examples).In a configuration of many particles, the coordinates of each interaction siteare specified by the particle’s fixed geometry (based on particle type), the positioncoordinates of the particle’s center of mass, and an orientation vector. The orienta-tion vector gives the particle’s orientation in the global coordinate system relativeto the orientation defined in the particle’s fixed geometry particle. The orienta-tions are represented using quaternions,119 which are vectors with four elementsused to describe an orientation with three degrees of freedom. The redundancyavoids singularities that would be encountered if Euler angles were used. With thebuilt in redundancy of the fourth element, quaternions must be normalized in or-der to uniquely specify an orientation. The orientation vector has four elements,Q = (q1, q2, q3, q4), where |Q| = 1.In the global configuration, each particle has a center of mass coordinateand orientation vector associated with it. A pair of rotation matrices are used toconvert between the particle coordinate system and the global coordinate system.The rotation matrix, A, transforms the fixed particle coordinates into the globalcoordinate system. The inverse of A transforms the global coordinates into the146local particle coordinates,119A =q21 + q22 − q23 − q24 2(q2q3 + q4q1) 2(q2q4 − q1q3)2(q2q3 − q1q4) q21 − q22 + q23 − q24 2(q3q4 + q2q1)2(q2q4 + q3q1) 2(q3q4 − q2q1) q21 − q22 − q23 + q24. (A.6)The coordinates of the interaction site α on ion i in the global configuration, riα,are calculated from the center of mass coordinates ri, the rotation matrix Ai, andthe reference geometry of the particle type asriα = ri +Aisα , (A.7)where sα is a column vector from the center of mass of particle i to the interactionsite α in the fixed particle coordinate frame.A.3 Interaction Potentials, Forces, PressureTo approximate the bulk properties of the salts without dramatically increas-ing the necessary computational resources, the simulation cell is repeated (tiled) inall directions. Applying periodic boundary conditions (PBCs) prevents surface in-teractions.119 The cell boundaries do not exert forces on the particles or disrupt themotions of the particles. When a particle moves through one side of the simulationcell, the particle coordinates are translated by the length of the cell so that it ap-pears to re-enter the simulation cell from the opposite side. A snapshot with PBCsis shown in Fig. 2.6 for a two-phase simulation.147Interaction PotentialsThe interaction potential and the definitions of terms within it was introducedin Chapter 2.2. The total potential energy is calculated for the configuration with theknown coordinates of each interaction site. The two types of interaction potentialsemployed here are the Lennard-Jones (LJ) potential and electrostatic potential. TheLJ potential is usually implemented in shifted and truncated form,uLJ(rij) =4εij[(σijrij)12−(σijrij)6]− uLJ(Rcut), if rij ≤ Rcut0, if rij > Rcut,(A.8)where rij = |rij| = |ri− rj|, and Rcut is set as a parameter. The entire LJ potentialis shifted up by uLJ(Rcut) and truncated at rij = Rcut. The shift eliminates thediscontinuity in the potential and any unphysical changes in the forces resulting froma particle moving across the Rcut threshold. For a collection of particles interactingthrough a LJ potential, the total potential energy is calculated asULJ =12N∑iN∑j 6=iuLJ(rij) (A.9)where the factor of 1/2 eliminates the double counting that occurs by summing overtwo particle indices i and j.The shifted and truncated LJ potential, as written in Eq. (A.8), neglects theinteractions at distances greater than Rcut. In order to get accurate energies andpressures in a simulation, a long-range correction term is added to the energy andpressure terms that accounts for the shifted and truncated potential. The correctionterm takes the form of an integral from Rcut to infinity, and assumes a uniform148density of LJ particles. The equation for the long-range correction of the potentialis given in Eq. (A.10), and the calculated quantity is added to ULJ from Eq. (A.9).The correction is calculated for the entire system of N particles in volume V ,U tail corrLJ =8piN23V εijσ3ij[13( σijRcut)9−( σijRcut)3]. (A.10)The long-range correction to the pressure is calculated asP tail corrLJ =16piN33V 2 εijσ3ij[23( σijRcut)9−( σijRcut)3], (A.11)which is added to each diagonal element of the stress tensor. In both the MD codewritten for this project and the Gromacs runs, long-range corrections are applied tothe potential energies and pressures.The electrostatic potential is long-ranged and truncating the interactionsbeyond the LJ potential cutoff distance, Rcut, would introduce spurious errors intothe potential and the force. The contributions to the potential energy are non-negligible and must be included in the calculation of the total potential and forces.The total electrostatic potential energy of point charges in a simulation cell withPBCs applied isUES =12N∑i,j=1′∑nqiqj|rij + nh|, (A.12)where n is a row vector of integers giving the translational directions of the im-age cells. In both programs, the maximum of each element of n is calculated asmax(nψ) = |Rcut/hψψ|, rounded up to the nearest integer, where ψ = 1, 2, or 3.The range for each element spans from −max(nψ) to max(nψ). The prime on then sum indicates that the i = j term of the n = (0, 0, 0) sum is omitted to prevent149calculating self-interactions. The prefactor of 1/2 is to remove the double countingthat occurs with the separate sums over i and j.The electrostatic potential is inversely proportional to the distance betweenpoint charges, and so the potential varies greatly at short distances and slowly atlarge distances. In order to account for these long-range, slowly varying interactions,the electrostatic potential is calculated using the Ewald sum, where the 1/r part ofthe potential is modified by dividing it into short- and long-range parts,119,136–1381r =f(r)r +1− f(r)r . (A.13)The modification, adding and subtracting a function f(r) to the potential, is de-signed to reduce the short-range contributions to negligible amounts at the trunca-tion limit (Rcut). A common, but by no means exclusive, choice for the functionf(r) is the complementary error function,erfc(x) = 2√pi∫ ∞xexp(−Ω2)dΩ, (A.14)where the error function iserf(x) = 1− erfc(x) = 2√pi∫ x0exp(−Ω2)dΩ. (A.15)The pair of functions is naturally substituted into the expression for 1/r where Ω isreplaced with κr, where κ is a parameter that set to ensure both parts of the sumconverge rapidly,1r =erfc(κr)r +erf(κr)r . (A.16)The first term on the right in Eq. (A.16) is the short-range part of the Ewald150summation. The function erfc(κr)/r decays exponentially as r → ∞ and has asingularity at r = 0. The short-range contribution is calculated in real space and istruncated for computational efficiency. The second term on the right of Eq. (A.16)has a long-range, slowly decaying tail. The long-range part is solved by takingthe Fourier transform of the modified function and summing over reciprocal spacevectors. Because of the sum over reciprocal space cell vectors, a self-correction termis applied that removes the contributions from a charge interacting with itself. Formolecules with multiple interaction sites, a molecular correction is applied to removethe intramolecular self-energy. The Ewald sum is conditionally convergent.120 Itis important to ensure that the simulation box is exactly neutral and that thecontributions from enough wave vectors are calculated in the reciprocal space sumto verify that the total energy has converged. The contributions to the Ewaldsum that take into consideration PBCs and non-cubic simulation cells are given inEq. (A.17) to Eq. (A.20).U real spaceES =18pi0N∑i,j=1′∑nqiqjerfc(κ|rij + nh|)|rij + nh|(A.17)U reciprocal spaceES =12piV′∑kexp (−pi2|r2k|/κ2)|r2k|· (A.18)N∑i=1∑α([ qiα√4pi0sin(2pirk · riα)]2+[ qiα√4pi0cos(2pirk · riα)]2)U selfES =−14pi0κ√pi∑i∑αq2iα (A.19)Umolecular correctionES =−18pi0κ√pi∑αqα∑βqβerf(κ|s|)|s| (A.20)In Eq. (A.18), k is a row vector of integers, spanning the reciprocal latticevectors in an analogous manner to n spanning real space lattice vectors. The vector151rk is calculated as rk = kh−1 where h−1 is the inverse of the cell matrix. Thevariables riα and qıα are the Cartesian coordinates and charge of site α on particlei.The number of k vectors retained must be sufficient for the reciprocal spacepart to converge, but keeping excess of k vectors increases computational cost. Theparameter κ controls the distribution of charge between the real space and reciprocalspace part. It is advantageous to use a value of κ that keeps the distribution even,to improve computational efficiency and to keep the relative error in both partsof the sum approximately equal. The parameter κ is usually optimized first bykeeping a large number of k vectors. Once κ is optimized, the number of k isreduced until reducing it further affects the convergence of the reciprocal spacesum (within a desired tolerance). The implementation of the Ewald sum can bechecked by computing the Madelung constants for different crystal lattices to adesired tolerance.Typically, another correction term is required to account for the surface ef-fects of building up the central simulation cell as a sphere in a medium with arelative permittivity, s. If that medium is a vacuum (s = 0), the correction mustbe applied to the potential.119 However, the appropriate boundary conditions for theionic systems we consider here are the so-called tinfoil boundary conditions, wheres →∞ and the correction term is zero.In Gromacs, the electrostatic interactions were calculated using the particle-mesh Ewald technique. The short-range interactions are calculated in the samemanner as the typical Ewald sum. The treatment of the long-range interactions isslightly different. The reciprocal space charge density is mapped on to a mesh witha certain grid spacing. The potential is calculated using a Fast Fourier Transform.152Forces and PressureThe forces on each particle are used to propagate MD simulations throughtime, the details are presented further below in Section B.1. The force is the negativegradient of the potential, f = −∇u, and is calculated in the evaluation routine ofboth the GA and MD programs. For the LJ pair potential, the force expressionis136–138fLJ(rij) = 24εij[2(σijrij)12−(σijrij)6]rij|r2ij|. (A.21)For the Ewald sum, the force contributions from the real space and the reciprocalspace terms aref real spacei =qi4pi0N∑j=1,j 6=iqj∑n[ 2κ√pi exp(κ|rij + nh|2)+erfc(κ|rij + nh|)|rij + nh|]rij + nh|rij + nh|2(A.22)f reciprocal spacei =18pi20V∑k 6=0exp(−pi2r2kκ2)·[qi sin(rk · ri)(N∑j=1qj cos(rk · rj))−qi cos(rk · ri)(N∑j=1qj sin(rk · rj))]rk|r2k|. (A.23)The self-term and the molecular correction in the Ewald sum will not contribute tothe force, as they are independent of position. The total force on particle i is thenused to advance its position in the MD program.For a particle with interaction site(s) off-center, the forces can result in atorque. It is convenient to advance the equations of motion in the particle’s local153coordinate frame (denoted with a superscript l) rather than the global coordinatesystem (denoted with a superscript g). The torque on a particle’s center in the localcoordinate frame, τ li , is calculated asτ li = Ai∑αsgα × f gα (A.24)where Ai is the rotation matrix from Eq. (A.6) to convert from the global to localparticle coordinates, sgα is the position vector from the ion center to site α in theglobal coordinate frame, and f gα is the force on site α in the global coordinate frame.The torque in the local frame is used in the equations of motion for the MD programdescribed in Appendix B.The virial theorem is used to calculate the instantaneous pressure tensor asPψν =1V(N∑i=1miviψviν +N∑i=1riψfiν), (A.25)where mi, vi, and fi are the mass, velocity, and force on particle i. The second termon the right in Eq. (A.25) with the product rifi must be calculated with care forcharged particles under periodic boundary conditions in non-cubic cells.136–139 Thehydrostatic (bulk) pressure is calculated as the average of the diagonal elements ofthe pressure tensor.154Appendix BMolecular DynamicsB.1 Time EvolutionMD simulations are propagated through time by integrating the equations ofmotion over small time steps. Two common methods of integrating time are the leap-frog algorithm119,120 and the Gear Predictor-Corrector scheme.119 The simulationsreported in this thesis were done using the leap-frog algorithm in Gromacs. In theleap-frog algorithm, the particle velocities and coordinates are advanced half of atime step out of sync,vi(t+ ∆t2)= vi(t− ∆t2)+ fi(t)mi∆t (B.1)ri (t+ ∆t) = ri(t) + vi(t+ ∆t2)∆ . (B.2)The particle velocities at time t, are calculated as the average of vi(t− ∆t2)andvi(t+ ∆t2).155A fifth-order Gear Predictor-Corrector scheme119 was used to propagatethrough time in the MD code produced for this project. The predictor-correctorscheme is described here using a stand-in variable Z. To propagate Z throughtime, Z, and its first four derivatives at time t, are advanced to predicted valuesat time t + ∆t using a Taylor series expansion. The superscript P on Z denotes apredicted value,ZP (t+ δt) = Z(t) +(dZdt)δt+ 12(d2Zdt2)δt2 + 16(d3Zdt3)∆t3 + 124(d4Zdt4)δt4.(B.3)In the Gear Predictor-Corrector scheme, the Taylor series expansion is also appliedto the time derivatives of Z. The set of prediction equations in matrix form isZP0 (t+ ∆t)ZP1 (t+ ∆t)ZP2 (t+ ∆t)ZP3 (t+ ∆t)ZP4 (t+ ∆t)=1 ∆t 12 ∆t2 16 ∆t3 124 ∆t40 1 ∆t 12 ∆t2 16 ∆t30 0 1 ∆t 12 ∆t20 0 0 1 ∆t0 0 0 0 1ZP0 (t)Z1(t)Z2(t)Z3(t)Z4(t)(B.4)where the subscript indicates the nth time derivative of Z. The predicted configura-tion (using ZP (t+ ∆t)) is evaluated and then corrected using equations of motion.In an example Newtonian system, the forces are evaluated at t+ ∆t using the posi-tions and interaction potential(s), and then the positions and their time derivativesare corrected usingr¨i = fi/mi. (B.5)For particles with more than one interaction site, the orientational equations of156motions are corrected in the local coordinate frame as119ω˙l,Cψ =(τ lψ + ωl,Pν ωl,Pφ (Iνν − Iφφ))/Iψψ, (B.6)where ω˙l,Cψ is the corrected time derivative of the angular velocity in local coordinatesand I is the principle moment of inertia matrix. The indices ψ, ν, and φ are axeslabels x, y, or z in the local coordinate system.The updated angular velocities are used to correct the orientation vector asq˙C1q˙C2q˙C3q˙C4= 12q1 −q2 −q3 −q4q2 q1 −q4 q3q3 q4 q1 −q2q4 −q3 q2 q10ωlxωlyωlz. (B.7)The difference between the predicted and corrected quantities, ∆Z, is thenused to correct the other time derivatives asZC0ZC1ZC2ZC3ZC4=ZP0ZP1ZP2ZP3ZP4+C0C1C2C3C4∆Z (B.8)where ∆Z = ZC −ZP is the difference between the predicted and calculated valueof Z from the evaluation routine, and the C matrix contains a set of numericalcoefficients (discussed further below). In the Newtonian example, ∆Z = f/m− a,157and applying the correction to the positions and their time derivatives is whatmaintains the Newtonian mechanics of the system.The matrix of numerical coefficients, C, are determined by the order of thedifferential equation of motion.119 The positions are governed by a second orderdifferential equation, whereas the angular velocities and orientations are governedby a first order differential equation. The values of the non-zero coefficients wereestablished by Gear.119Table B.1: Numerical coefficients used in the Gear Predictor-Corrector algorithm.119Coefficient 1st order DE 2nd order DEC0 251/720 19/120 or 19/90C1 1 3/4C2 11/12 1C3 1/3 1/2C4 1/24 1/12The Gear predictor-corrector algorithm was used to advance the positions,orientations, Nose´-Hoover thermostat, and the cell dynamics in the MD programwritten for this project.B.2 Simulation ConstraintsB.2.1 Temperature — Particle MotionA Nose´-Hoover thermostat100–102 regulates the temperature of an ensembleof particles by adding a correction term to the equations of motion. The particle158accelerations are propagated through time asai =fimi− ξvi, (B.9)where ξ is a friction coefficient that couples the ensemble of particles to the desiredtemperature, Tset. The friction coefficient, ξ, propagates through time asdξdt =1Q[N∑imi |vi|2 − 3NkBTset], (B.10)where kB is the Boltzmann constant. The term in brackets measures the differencebetween the instantaneous temperature in the simulation and the target temper-ature. The temperature is regulated through ξ, which increases or decreases theparticle motion when there is a difference between the desired and actual temper-ature. A similar modification is made to Eq. (B.6) to thermostat the rotationalmotion, where a separate friction coefficient times the angular velocity is subtractedfrom the right hand side.The response time of the thermostat is determined by the “heat bath mass,”Q, which is a parameter that must be chosen with care. A small Q gives a fasterresponse time and has a larger influence on the particle dynamics. A loose thermalcoupling (large Q) will have a slow response time, meaning that on one hand, anout-of-equilibrium system may take a long time to reach equilibrium but, on theother hand, the perturbations of the particle dynamics are not as large. For theMD code produced for this project, Q = 5 is an appropriate value. In Gromacs, thestrength of the thermal coupling for all thermostats is set in terms of a relaxationtime, τT . For the Nose´-Hoover thermostat, τT ∝ Q1/2. In the Gromacs simulationsdone here, the relaxation time was 0.1 ps, with separate coupling groups for the159cations and anions.B.2.2 Pressure — Simulation CellThe simulation cell size and shape are both allowed to fluctuate in our NPTmolecular dynamics simulations. The desired pressure, Pset was maintained usingthe Parrinello-Rahman barostat103,104 in the Gromacs simulations. The equationsof motion use fractional coordinates, rFi = (h−1)Tri,r¨Fi =(h−1)Tfimi− g−1g˙r˙Fi (B.11)h¨T = 1W (P − IPset)V h−1, (B.12)where the single and double dots over variables represent the first- and second-derivatives with respect to time, and I represents the 3 × 3 identity matrix. Here,we see the dependencies on the simulation cell vectors (h), the metric tensor g, andthe instantaneous pressure tensor, P . The effective cell mass, W , determines theresponsiveness of the cell dynamics. A slow response (large W ) means the size andshape of the container will fluctuate slowly, whereas, a small W will lead to a moredynamic container. The relaxation parameter W of the barostat is analogous to therelaxation parameter Q of the thermostat.In the code produced for this project, the size and shape fluctuations wereimplemented using Souza and Martin’s metric tensor approach.139 In this approach,the equations of motion for the particle dynamics are the same as Eq. (B.11). Theparticle dynamics are governed by the same equation for both the metric tensorbarostat and Parrinello-Rahman barostat, but trajectories from identical initial con-figurations would produce different particle motions. The dynamics of the cell affect160the particle trajectories through g˙ in Eq. (B.11), and in the metric tensor approach,the cell dynamics are not governed by Eq. (B.12), but byg¨γδ =12W√det g( Pγδ√det g − Psetgγδ)+(g˙γζgζηg˙ηδ − gζηg˙ζηg˙γδ)+ 12(g˙ζηgηλg˙λνgνζ)gγδ, (B.13)where the Einstein notation has been used to signify summation over repeated Greekindices. Details for the implementation can be found in the paper by Souza andMartins.139161Appendix CGenetic AlgorithmGenetic algorithms optimize multi-parameter problems using concepts bor-rowed from evolutionary biology such as survival of the fittest.140 The algorithmstarts by taking a collection of possible solutions, called a population, and repeats acycle of evaluation and evolution until the conditions for termination are satisfied.In each cycle, the possible solutions are ranked by set criteria and the algorithm isdesigned to optimize the candidate structures.Genetic algorithms have been developed and used for finding crystal struc-tures and cluster geometries.114,141–148 In the present implementation, a candidatestructure is a viable crystal unit cell and the GA is designed to minimize the poten-tial energy.An overview of the GA program is given as part of Fig. A.1. The generationof the initial set of candidate crystal structures, their evolution through applyinga variety of mutation operators, and how the GA terminated are discussed in this162section.C.1 InitializationEach candidate crystal structure is specified by three unit cell vectors, a, b,c, three angles between those vectors, α, β, γ, and the number, type, fractionalcoordinates, and orientation (if applicable) of each particle within the unit cell. Theparticle “type” specifies the particle label, the rigid geometry of the interactionsites, and the interaction parameters for each site within the particle. The numberof particles and the ratio of particle types is set before the run to reduce the numberof variables the GA has to optimize. To sample different unit cell possibilities, thenumber of particles within the unit cell is varied over separate runs of the GA.The GA is used for two types of systems here, salts and binary LJ mixtures. Forthe salts, charge neutrality limits the possible particle type ratios, whereas for thebinary LJ mixtures, the ratio of particle types is varied as a run input.The six unit cell parameters are initialized as cubes with the length of eachside set to V 1/3i where Vi is an approximate desired volume of the unit cell (near areduced density of 1.0). An option was included to keep all cell angles fixed at 90◦for the entire GA run to exclusively sample orthorhombic configurations.The initial particle labels (types) are assigned with the help of a pseudo-random number (PRN) generator. The range [0,1] is divided into the number ofparticle types present. A PRN on [0,1] is generated for each particle (in serial), andthe particle label is assigned based on comparing the label ranges and the PRN.Counts of the assigned labels are kept to ensure the assigned particle types matchedthe input specifications. The fractional coordinates along each axes are assigned163using the PRN on [0,1] for each particle. Each particle is initialized with the sameorientation.C.2 Configuration ChecksBefore each evaluation of the energies of the population of crystal structurecandidates, the configurations are put through some “rational” checks. The candi-dates are checked for obvious signs that the crystal structure will yield a poor (high)energy or unphysical structure.The first check screens the unit cells for lengths that are too short or toolong and angles that are too acute or too obtuse. The cell length limits are com-puted from the particle type with the largest LJ diameter, σmax, and the numberof particles within the cell. The minimum cell length allowed is Nσmax/6 and themaximum length allowed is Nσmax/2. If the cell length falls outside of these bounds,a replacement cell length is generated using a PRN on [Nσmax/6,Nσmax/2]. Eachangle has to be between 45◦ and 120◦, otherwise the cell is reset.The fractional coordinates of each particle must be between 0.0 and 1.0, andare translated by ±1.0 as necessary. The fractional coordinates are also checked foroverlaps. If the distance between two particles along all three axes are all below athreshold of 0.05 (in fractional units), then one particle’s fractional coordinates areregenerated using the PRN generator, and the configuration is checked again.After the entire population passes the checks, the energies are evaluated(as discussed in Appendix A.3). Before the mutation operators are applied, thecandidate structures are ranked from best to worst by their potential energies.164C.3 Configuration EvolutionAside from the elitism operator, which preserves the previous generationsbest candidates, each mutation operator is designed to change some attribute ofthe candidate crystal structure. A collection of 15 mutation operators that balanceboth drastic changes and minor tweaks have been implemented in order to effectivelysearch for the best crystal structure.The rates of application of each mutation operator are set as program param-eters, and have been tested to ensure a robust collection of candidates are producedand that the GA is able to find, and converge to, the same crystal structure frommultiple runs with different initializations of the PRN generator. The mutationoperators are applied to candidates selected using the PRN with the best candidateexcluded.ElitismThe elitism operator creates a perfect copy of the best candidate crystalstructures from the current generation and places them into the next generation.The sole function of this operator is to preserve the best candidates through gen-erations. For a population of 50 candidate structures, the 6 best structures arepreserved using the elitism operator (12%). Keeping the elitism operator near 10%maintains the diversity of the population long enough to explore a variety of crystalstructures while still converging within a reasonable number of generations.165C.3.1 Diversity OperatorsThe operators described in this section are particularly successful at produc-ing lower energy configurations in the early stages of a run, because the new config-urations (offspring) are very different from the initial configurations (parents). Eachof these operators is used to create 5 — 10% of the next generation.CrossoverThe crossover operator splices two candidate structures together to createtwo “offspring.” The attributes of the first parent are copied to the first offspringuntil the crossover point is reached. After the crossover point, the attributes of thesecond parent are copied to the first offspring to complete the attribute set of thefirst offspring. The process is repeated to create a second offspring where the orderof the parents is switched (the crossover point remains the same).The cell lengths and angles are preserved from the dominant parent. Thecrossover point is determined by a PRN generated on [2,N-1]. The recombinationis applied to the particles only. The particle labels, fractional coordinates, and theparticle orientations are transferred from one parent before the crossover point andthe other parent after.BlendThe blend operator, like the crossover operator, creates two new candidatestructures from two parents in the current generation. In this operator, the frac-tional coordinates of the two parent configurations are sorted and then paired (onefrom each parent) by closest distance. The offspring’s fractional coordinates are then166generated by taking a weighted average of the positions. The average is weighted2/3 to the dominant parent and 1/3 to the recessive parent. Switching the domi-nant/recessive roles of the parents produces a second offspring configuration. Theparticle labels and orientations in the offspring configurations are directly copiedfrom the dominant parent. The cell lengths and angles are also subject to theweighted average.Geometric MeanThe geometric mean operator produces one offspring from two parents. Eachof the cell parameters and fractional coordinates are generated by combining thecorresponding attributes of both parents in a geometric fashion. For example, thelength of cell vector a is computed as |achild| =√|aMom| |aDad|. The particle labelsand orientations are all reassigned using the PRN generator.SwapThe swap operator mutates a configuration by switching two fractional coor-dinates of two different particles. The two particles and which fractional coordinatesare to be swapped are selected using the PRN. An example would be swapping thefractional coordinate along b on particle 3 with the fractional coordinate along a ofparticle 5.RippleThis operator is constructed from the description of the Ripple operator fromthe excellent work of D.C. Lonie and E. Zurek.148 The ripple operator displaces the167coordinates in one direction based on the position of the other two coordinates andhow many wavelengths are chosen to fit in the cell. The maximum wavelength isset to 2, and the magnitude of the ripple is set to 0.15 in this work. This operatorcopies the cell attributes, particle labels, and orientations directly from the parentconfiguration. Only the fractional coordinates of the particles change.StrainThe strain operator also follows from the work of D.C. Lonie and E. Zurek.148The fractional coordinates, orientations, and labels of the particles of the parentconfiguration are copied exactly to the offspring. The parent’s cell parameters areadjusted by applying a symmetric strain matrix that includes six PRNs that skew(or deskew) the cell vectors.Swap CellThe swap cell operator generates two offspring from two parents. The prop-erties of the particles are transmitted from one parent to one offspring, but theunit cell parameters are taken from the other parent. There are no mutations ofthe attributes, but the two new combinations of attributes create offspring that aredistinct from the parents (unless, of course, the cell parameters of the two parentsare similar).Replace AttributesTwo operators that simply replace part of the parent configuration with a newone that is generated pseudo-randomly. One operator replaces the three fractional168coordinates and orientations of a single particle in the parent configuration. Anotheroperator replaces the orientations of all particles within a configuration. Thesereplacement operators are applied least of the diversity operator set—they eachcreate about 4% of the next generation.C.3.2 Refinement OperatorsThe changes introduced by the refinement operators are typically small.Whereas the diversity operators discussed earlier are more successful at the be-ginning of a GA run, these refinement operators have more success at improvingthe configurations near the final stages of the GA run. Each of the following opera-tors are used to generate less than 5% of the configurations for the next generationexcept for the adjust volume operator (which is closer to 10%).Adjust VolumeThe adjust volume operator scales the cell lengths by a maximum of 5%.Two PRNs are used in this operator. One determines the magnitude of the vol-ume shift (generated on [0,1] and multiplied by the maximum volume adjustmentpermitted), and the second PRN (generated on [-1,1]) determines the sign of thevolume fluctuation.ReorientThe reorient operator changes the orientations of all of the particles from theparent configuration. The orientations are rotated by a maximum of 10◦ about oneaxis. Each particle is rotated by a different amount. The parent’s cell parameters,169fractional coordinates, and particle labels are preserved.Adjust Cell AngleThe adjust cell angle operator changes one cell angle by a maximum of 20◦.The actual angle and the direction are generated using PRNs in a manner similarto the volume adjustment operator. The cell angle that is adjusted (α, β, or γ) isselected pseudo-randomly.Adjust Cell LengthsThree operators, in addition to the volume adjust operator described above,are used to manipulate the cell axes lengths. The first operator makes all threecell lengths equal by taking an average of the parent’s cell lengths, which is usefulfor accessing cubic structures. The second operator scales a single axis (selectedpseudo-randomly) by a maximum of 10%, in a manner similar to the volume adjustoperator. The third operator contracts the longest cell axis by up to 15%.C.4 TerminationThe GA terminates after executing a predetermined number of generations(usually 200), regardless of how long the “best” crystal structure has been ranked#1. Another common termination scheme stops the run once the energy of the bestcrystal structure has not been improved upon for a set number of generations.In order to be confident that the final configuration is one of the lowest energycrystal structures, the GA is run multiple times with different PRN initializations170and different parameters (particularly N , the number of particles in the unit cell).C.5 GA ResultsThe program was tested on binary mixtures of LJ particles. The two particletypes, A and B, have LJ length parameters σAA = 1.0, σAB = 0.8, and σBB = 0.88,while the LJ energy parameters are εAA = 1.0, εAB = 1.5, and εBB = 0.5. Theseparameters do not follow the Lorentz-Berthelot mixing rules119 and are designed tofavour the mixing of the two particle types. The parameters roughly correspond tothe nickel phosphorus (Ni−P) system.149,150 The A80B20 composition correspondsto a eutectic in the Ni−P system, which has not been reproduced in the LJ binarymixture.151 The GA results (the structures, densities, and potential energies atT = 0 and P = 0) agree with those found in the literature.151The GA was not used for the model salts studied here because the majority(about 84%) of the salts spontaneously crystallized upon cooling in MD simulations.An alternate approach to find crystal structure candidates (preparing and heatingcrystal structures in MD simulations) was used because it was convenient for thesimple salts studied here.171Appendix DTemperature TablesTable D.1: The five characteristic melting-related temperatures for the size sym-metric 2L - 1C salts that crystallize spontaneously. The Tm values are the esti-mates from the hysteresis loops [Eq. (2.5)]. The two-phase simulations give thebounding temperatures for the CsCl-solid to liquid transition. All temperaturesare in kelvin.Cation Hysteresis Loops Two-phaseT− T+ Tm Ts Tl1C 1050 1400 1238 1250 13002L33-06 1050 1400 1238 1250 13002L33-08 1050 1400 1238 1250 13002L33-10 1050 1300 1182 1200 12502L33-12 1050 1300 1182 1200 12502L33-14 1050 1250 1154 1150 12002L33-16 950 1300 1139 1150 12002L33-18 1000 1150 1078 1100 11502L33-20 900 1150 1032 1050 11002L33-22 850 1200 1040 1050 11002L33-24 800 1000 906 950 1000Continued on next page. . .172Table D.1: (Continued)Cation Hysteresis Loops Two-phaseT− T+ Tm Ts Tl2L50-06 1050 1250 1154 1250 13002L50-08 1000 1350 1188 1200 12502L50-10 1050 1250 1154 1200 12502L50-12 950 1300 1139 1150 12002L50-14 950 1350 1167 1100 11502L50-16 850 1150 1011 1050 11002L50-18 750 1000 884 950 10002L50-20 700 950 835 850 9002L50-22 700 850 778 750 8002L50-24 650 800 729 750 8002L67-06 1000 1450 1246 1200 12502L67-08 1000 1350 1188 1200 12502L67-10 950 1300 1139 1150 12002L67-12 850 1150 1011 1100 11502L67-14 750 1100 942 1000 10502L67-16 700 900 806 850 9002L67-18 550 850 716 800 8502L67-20 550 850 716 750 8002L100-06 950 1250 1110 1150 12002L100-08 900 1200 1060 1100 11502L100-10 700 1250 1015 1000 1050173Table D.2: The five characteristic melting-related temperatures for the 2L - 1Csalts that vitrify. Here, T− is the glass transition temperature, Tg, as definedearlier, and T+ is for the melting of the 111n crystal structure. The temperaturesfrom the two-phase simulations are for the 111n solid to liquid transition. Alltemperatures are in kelvin.Cation Hysteresis Loops Two-phaseT− T+ Tm Ts Tl2L67-20 400 850 667 700 7502L67-22 400 850 667 650 7002L67-24 400 800 634 600 6502L100-12 450 900 714 700 7502L100-14 400 1000 768 750 8002L100-16 400 1000 768 800 8502L100-18 400 1000 768 800 8502L100-20 400 950 734 600 6502L100-22 400 800 6342L100-24 350 650 523174Appendix ECrystal StructuresE.1 111n StructureThe 111n structure was shown and described in Chapter 3 (Fig. 3.3).Space Group 111 D2d-1 P-42mOrigin at 0.23450 0.24000 0.82550Values of a,b,c,alpha,beta,gamma:0.9420 0.9420 1.0123 90.000 90.000 90.000Atomic positions in terms of a,b,c:CM1 Wyckoff position n, x = -0.23513, z = 0.1775CC1 Wyckoff position n, x = -0.23538, z = 0.3390AM1 Wyckoff position n, x = -0.23550, z = -0.3725Figure E.1: Selected output from running the Findsym program on the unit cellabove. The tolerance was set to 0.03.175Table E.1: The contents of a 111n unit cell from simulations of the 2L100-16 - 1Csalt. The table gives the unit cell dimensions, site labels, and xyz coordinates. Thecoordinates and the unit cell dimensions are in nanometers.unit cell lengths (a,b,c) 0.93929 0.94383 1.01230unit cell angles (α, β, γ) 90.0 90.0 90.0Site label x y zCM1 0.000 0.480 0.650CC1 0.000 0.477 0.490CM2 0.469 0.472 0.001CC2 0.470 0.474 0.170CM3 0.000 0.000 0.000CC3 0.000 0.000 0.157CM4 0.469 0.009 0.641CC4 0.469 0.007 0.481AM1 0.000 0.473 0.202AM2 0.470 0.480 0.457AM3 0.000 0.003 0.448AM4 0.469 0.005 0.193176E.2 Var-CsCl StructureThe var-CsCl crystal structure was found as a result of the earlier work onthe A100 2L67-24 - 1C salt, but was not fully investigated then. The structure isshown in Fig. E.2, and is an orientationally ordered variation on the CsCl crystalstructure.(a) yz-view(b) xy-view (c) xz-viewFigure E.2: Orientationally ordered crystal structure of the A100 2L67-22 - 1C salt.The cation centers, cation off-center sites, and anions are shown in blue, white, andred, respectively, and are not shown to scale.177E.3 Var-111n StructureThe var-111n structure is a variation of the 111n crystal structure. The crystalstructure adjusts to accommodate cations that are slightly larger than the anions.The B133 2L67-22 - 1C and B133 2L67-26 - 1C salts are the two that rearrangedinto the var-111n structure.Figure E.3: The var-111n crystal structure shown for the B133 2L67-22 - 1C salt.The cation centers, cation off-center sites, and anions are shown in blue, white, andred, respectively, and are not shown to scale.178Appendix FInput FilesExample input files for heating the 111n crystal of the A100 2L100-16 -1C salt are included in this Appendix. The execution script (bash) is a *.pbs filedesigned for job submission on Orcinus and the example included here will heat aninitial configuration from an initial temperature of 50 K up to 1200 K (or until thewall clock limit is reached, whichever is first). The *.pbs files were generated in sets(cooling, heating, and two-phase simulations) for each salt using a bash script.The *.gro file is the initial configuration file that gives the particle labels,positions, and optionally, velocities, and the lengths that define the simulation cell.The initial *.gro files were crystal structures (either CsCl at a reduced density of0.8 or a variety of crystal structures at a reduced density of 1.1, depending on thedesired target phase) and were generated using a small program I wrote.The *.top file specifies the topology and number of each type of particle inthe system. (Note that the number of each particle type specified in the *.top file179and given in the *.gro file must match.) The *.top files vary for each salt and weregenerated using a bash script and were automatically modified before the simulationto ensure the numbers of particles matched the *.gro input file.The *.mdp file specifies the run parameters for Gromacs. This particular*.mdp example file is for a 4.0 ns NPT run (with the barostat parameters used in thehysteresis loops) at a temperature of 600 K. I had folders of *.mdp files symbolicallylinked in each running folder for simulation type. The hysteresis loops required setsof *.mdp files for each cutoff distance. The two phase simulations required foursets of *.mdp files for each cutoff distance (energy minimization, NV T , anisotropicNPT , and isotropic NPT ). The isotropic NPT .mdp files were not the same as theisotropic NPT .mdp files for the hysteresis loops; the barostat relaxation parameterwas slower for the two-phase simulations. Within each folder, there was an *.mdpfile for each simulation temperature.180#!/bin/bash#PBS -S /bin/bash#PBS -l procs=8#PBS -l walltime=80:00:00#PBS -r n#PBS -N A100_2L100-16_1C_H___111n#PBS -M <<email address>>#PBS -m eacd $PBS_O_WORKDIRecho "Current working directory is `pwd`" module load gromacs/4.6.2# *******# ion specifics# *******# cationgeom=2Lgeom1=100cdi=16cda=$(printf %02d $cdi)# aniongeom2=1C# ion sizesSize_ratio=100Cat_Size=050Ani_Size=050# *******# run specifics# *******rcutoff=170diri="H"# diri is a letter to indicate the temperature direction, H = heat a prepared crystal structure from 50 K.folder1= #source folder for .mdp foldersfolder2= #destination folder for runs# labelssalt=$geom$geom1"-"$geom2label=$geom"-"$cda"_"$geom2label2=$geom$geom1"-"$cda"_"$geom2longlabel=$geom$geom1"-"$cda"_"${Cat_Size}"_"$geom2"_"${Ani_Size}ln -s <<path to $folder1>>/npt_iso/npt_$rcutoff npt temp1=1200temp2=50grofile="50.gro"# "50.gro" is the name of the input .gro file for the run at the first temperature. For the heating runs (diri="H"), it is a perfect crystal structure set at a reduced density of 1.1.Figure F.1: Example .pbs file (page 1/3)181# *******# Main temperature loop:# *******((temp = $temp2)) while ((temp<=temp1)) do   if [ -s $grofile ]  then    ndxfile=$label2"_"$temp".ndx"    echo -e "aCM\naCC\naAM\nq\n" | make_ndx_mpi_d -f $grofile -o $ndxfile    nptfile="npt/npt"$temp".mdp"    topfile=$longlabel".top"    tprfile=$temp".tpr"    # make sure .top and .gro molecule numbers match. Trust the .gro file.    declare -i ntot=$(sed -n 2p $grofile)    nmol=$(( ntot * 2 / 3 ))    ncat=$(( ntot / 3 ))    nani=$(( ntot / 3 ))    localtop="local_top.top"    head -n -2 $topfile > $localtop    echo " CAT        "$ncat >> $localtop    echo " ANI        "$nani >> $localtop    mv $localtop $topfile    # simulate system    grompp_mpi_d -f $nptfile -c $grofile -p $topfile -o $tprfile    mpirun -np 8 mdrun_mpi_d -s $tprfile    if [ -s traj.trr ]    then      # rename output files to contain the run temperature      mv confout.gro $temp"_out.gro"      mv ener.edr $temp"_ener.edr"      mv traj.trr $temp"_traj.trr"      mv mdout.mdp $temp"_mdout.mdp"      mv md.log $temp"md.log"      mv state.cpt $temp"state.cpt"      # analyze the trajectory       echo -e "1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 0" | g_energy_mpi_d -f $temp"_ener.edr" -nmol $nmol -b 2000      mv energy.xvg $temp"_energy.xvg"Figure F.1 Continued: Example .pbs File (page 2/3)182      GRO=$temp"_out.gro"      TRAJ=$temp"_traj.trr"      RDF1=$temp"_rdf1_CM_CM.xvg"      RDF2=$temp"_rdf2_CC_CC.xvg"      RDF3=$temp"_rdf3_CM_AM.xvg"      RDF5=$temp"_rdf5_AM_AM.xvg"      RDF7=$temp"_rdf7_CC_AM.xvg"      RDF8=$temp"_rdf8_CC_CM.xvg"      echo -e "4\n4\n" | g_rdf_mpi_d -f $TRAJ -s $GRO -n $ndxfile -o $RDF1 -b 3000      echo -e "5\n5\n" | g_rdf_mpi_d -f $TRAJ -s $GRO -n $ndxfile -o $RDF2 -b 3000 -cut 0.3      echo -e "4\n6\n" | g_rdf_mpi_d -f $TRAJ -s $GRO -n $ndxfile -o $RDF3 -b 3000      echo -e "6\n6\n" | g_rdf_mpi_d -f $TRAJ -s $GRO -n $ndxfile -o $RDF5 -b 3000      echo -e "5\n6\n" | g_rdf_mpi_d -f $TRAJ -s $GRO -n $ndxfile -o $RDF7 -b 3000 -cut 0.15      echo -e "5\n4\n" | g_rdf_mpi_d -f $TRAJ -s $GRO -n $ndxfile -o $RDF8 -b 3000 -cut 0.15      echo -e "4" | g_msd_mpi_d -f $TRAJ -s $tprfile -n $ndxfile -ten -o $temp'_CM_msd.xvg'      echo -e "6" | g_msd_mpi_d -f $TRAJ -s $tprfile -n $ndxfile -ten -o $temp'_AM_msd.xvg'      echo -e "4" | g_velacc_mpi_d -f $TRAJ -s $tprfile -n $ndxfile -o $temp"_"$diri"_CM_VACF.xvg" -b 2000      echo -e "5" | g_velacc_mpi_d -f $TRAJ -s $tprfile -n $ndxfile -o $temp"_"$diri"_CC_VACF.xvg" -b 2000      echo -e "6" | g_velacc_mpi_d -f $TRAJ -s $tprfile -n $ndxfile -o $temp"_"$diri"_AM_VACF.xvg" -b 2000      echo -e "2" | g_rotacf_mpi_d -f $TRAJ -s $tprfile -n $ndxfile -o $temp"_"$diri"_CAT_RotACF1.xvg" -P 1 -d -b 2000      echo -e "2" | g_rotacf_mpi_d -f $TRAJ -s $tprfile -n $ndxfile -o $temp"_"$diri"_CAT_RotACF2.xvg" -P 2 -d -b 2000      echo -e "3" | g_rotacf_mpi_d -f $TRAJ -s $tprfile -n $ndxfile -o $temp"_"$diri"_ANI_RotACF1.xvg" -P 1 -d -b 2000      echo -e "3" | g_rotacf_mpi_d -f $TRAJ -s $tprfile -n $ndxfile -o $temp"_"$diri"_ANI_RotACF2.xvg" -P 2 -d -b 2000    fi    grofile=$temp"_out.gro"  fi# change the temp (note that grofile is set to the final configuration of the run at the last temperature, so it is used as the initial configuration at the new temperature.)((temp+=50))done   # *******# tidy the run folder# *******rm -f step*mkdir gromkdir rdfmkdir energiesmkdir ACFmkdir msdmkdir runstuffmv *_out.gro gro/.mv *rdf*.xvg rdf/.mv *ener.edr energies/.mv *energy.xvg energies/.mv *VACF* ACF/.mv *RotACF* ACF/.mv *_msd.xvg msd/.mv *.log runstuff/.mv *.cpt runstuff/.mv *.mdp runstuff/.mv *.tpr runstuff/.mv *.ndx runstuff/.Figure F.1 Continued: Example .pbs file (page 3/3)1832L100-16_1C 2592    1CAT     CM    1   0.090   0.017   5.606    1CAT     CC    2  -0.070   0.008   5.603    2CAT     CM    3   0.448   0.489   0.008    2CAT     CC    4   0.608   0.490   0.007    3CAT     CM    5   0.452   0.018   0.495    3CAT     CC    6   0.612   0.017   0.493    4CAT     CM    7   0.101   0.455   0.427    4CAT     CC    8  -0.059   0.458   0.427    5CAT     CM    9   0.111   0.014   0.917    5CAT     CC   10  -0.049   0.019   0.916    6CAT     CM   11   0.478   0.497   0.965    6CAT     CC   12   0.638   0.498   0.961    7CAT     CM   13   0.481   0.030   1.439    7CAT     CC   14   0.641   0.028   1.436    8CAT     CM   15   0.129   0.519   1.401    8CAT     CC   16  -0.031   0.521   1.405    9CAT     CM   17   0.128   0.070   1.879    9CAT     CC   18  -0.032   0.071   1.879   10CAT     CM   19   0.485   0.490   1.903   10CAT     CC   20   0.645   0.489   1.906   11CAT     CM   21   0.492   0.017   2.380   11CAT     CC   22   0.652   0.012   2.380   12CAT     CM   23   0.132   0.531   2.346   12CAT     CC   24  -0.028   0.531   2.349...  860CAT     CM 1719   5.077   5.247   4.213  860CAT     CC 1720   4.917   5.248   4.215  861CAT     CM 1721   5.082   4.766   4.717  861CAT     CC 1722   4.922   4.765   4.716  862CAT     CM 1723   5.420   5.230   4.691  862CAT     CC 1724   5.580   5.229   4.692  863CAT     CM 1725   5.440   4.716   5.193  863CAT     CC 1726   5.599   4.719   5.191  864CAT     CM 1727   5.083   5.216   5.167  864CAT     CC 1728   4.923   5.215   5.171  865ANI     AM 1729   0.650   0.497   0.484  866ANI     AM 1730   0.625   0.004   0.013  867ANI     AM 1731   5.866   0.514   5.619  868ANI     AM 1732   5.879   5.649   0.432  869ANI     AM 1733   0.681   0.494   1.436  870ANI     AM 1734   0.661   0.022   0.964  871ANI     AM 1735   5.889   0.482   0.909  872ANI     AM 1736   5.896   0.056   1.398... 1717ANI     AM 2581   5.642   5.259   3.261 1718ANI     AM 2582   5.644   4.797   2.825 1719ANI     AM 2583   4.909   5.251   2.784 1720ANI     AM 2584   4.896   4.761   3.311 1721ANI     AM 2585   5.615   5.224   4.214 1722ANI     AM 2586   5.634   4.799   3.761 1723ANI     AM 2587   4.897   5.259   3.731 1724ANI     AM 2588   4.889   4.779   4.248 1725ANI     AM 2589   5.618   5.190   5.161 1726ANI     AM 2590   5.616   4.755   4.710 1727ANI     AM 2591   4.885   5.234   4.690 1728ANI     AM 2592   4.897   4.747   5.192   5.96489   5.66679   5.65433Figure F.2: Truncated Example of a .gro file. Many particles (Cations 13-859 andAnions 873-1716) have been omitted for brevity.184[ defaults ]; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ1 2 yes 1.0  1.0[ atomtypes ];name  atnum  mass    charge   ptype   sigma    epsilonCM     11    115.000  0.000      A    0.500000  3.613285CC      1      5.000  0.000      A    0.000000  0.000000AM     21    120.000  0.000      A    0.500000  3.613285[ nonbond_params ];i     j  func    sigma   epsilonCM     CM   1   0.500000  3.613285CM     AM   1   0.500000  3.613285AM     AM   1   0.500000  3.613285CM     CC   1   0.000000  0.000000CC     CC   1   0.000000  0.000000AM     CC   1   0.000000  0.000000[ moleculetype ];molname  nrexclCAT   2    [ atoms ];id at type res nr residu name at name cg nr charge 1       CM    1          CAT      CM   2      0.0000000000000002       CC    1          CAT      CC   2      1.000000000000000[ constraints ]; ai   aj funct        b0   1     2     2     0.16000     [ exclusions ]    1 22 1[ moleculetype ]; molname  nrexclANI   2[ atoms ]; id at type res nr residu name at name cg nr charge1       AM    1          ANI      AM   2 -1.000000[system ]2L100-16_050_1C_050[ molecules ] CAT        984 ANI        984Figure F.3: Example .top file185; Run parametersintegrator = md   ; leap-frog integratornsteps   = 4000000 ; 4.0 nsdt       = 0.001   ; 1 fscomm_mode   = linear  ; remove center of mass motionnstcomm     = 10comm_grps   = system; Output controlnstxout = 1000 ; save coordinates every psnstvout = 1000 ; save velocities every psnstenergy = 1000 ; save energies every psnstlog = 1000 ; update log file every ps; Bond parameterscontinuation = no ; first dynamics runconstraint_algorithm = lincs ; holonomic constraints constraints = all-bonds ; all bonds (even heavy atom-H bonds) constrainedlincs_iter = 1 ; accuracy of LINCSlincs_order = 6 ; also related to accuracylincs_warnangle = 90.0 ; much larger than normal.; Neighborsearchingns_type = grid ; search neighboring grid cellsnstlist = 5   ; update listrlist = 1.7 ; short-range neighborlist cutoff (in nm)rcoulomb = 1.7 ; short-range electrostatic cutoff (in nm)rvdw   = 1.7 ; short-range van der Waals cutoff (in nm); Electrostaticscoulombtype    = PME  ; Particle Mesh Ewald for long-range electrostaticspme_order    = 4    ; interpolationfourierspacing = 0.2  ; grid spacing for FFT; Temperature coupling is ontcoupl = nose-hoover   ; thermostat typetc-grps = CAT ANI     ; two coupling groups - more accuratetau_t = 0.1 0.1       ; time constant, in psref_t = 600 600     ; reference temperature, one for each group, in K; Pressure coupling is onpcoupl   = Parrinello-Rahmanpcoupltype  = isotropictau_p       = 5.0compressibility = 3e-5ref_p           = 1.0; Periodic boundary conditionspbc = xyz ; 3-D PBC; Dispersion correctionDispCorr = EnerPres ; account for cut-off vdW scheme; Velocity generationgen_vel = no ; assign velocities from Maxwell distributiongen_temp = 0 ; temperature for Maxwell distributiongen_seed = -1 ; generate a random seedFigure F.4: Example .mdp file186


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items