Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

A molecular dynamics investigation of the dissolution of molecular solids Anand, Abhinav 2017

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

Item Metadata


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

Full Text

A molecular dynamics investigation ofthe dissolution of molecular solidsbyAbhinav AnandB.Tech., I.I.T. Guwahati, 2015A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFMASTER OF SCIENCEinThe Faculty of Graduate and Postdoctoral Studies(Chemistry)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)September 2017c© Abhinav Anand 2017AbstractThe dissolution of molecular solids is an important process, which has been studied for overa century. However, a lot of work is still needed for a detailed understanding of the molecularmechanism of dissolution, because of the complex nature of many molecular solids, and thelarge time scales required for simulation studies. In this thesis we study the dissolution ofmolecular solids, to examine if classical models (which assume that the rate is proportionalto an active surface area) can be used to describe the dissolution profile of these solids.Urea and aspirin molecules are used as models, to study the dissolution process in waterunder sink conditions, because of their contrasting solubility in water. The dissolution rate indifferent water models was examined and it was found that they differ considerably. However,the overall mechanism for the dissolution process remains the same. Dissolution was foundto be an activated process with the detachment of molecules from the crystal being the ratelimiting step. Crystals with different shapes (cubic and cylindrical) were used to study theeffect of shape on the dissolution process.The dissolution process for urea was found to occur in three steps, an initial rapid stage,where the molecules at the edges and corners go into the solution, a long intermediate stagewith a nearly constant dissolution rate, and a final stage where the crystals lose their crys-talline structure and dissolve completely. The fixed rate law stage was found to be describedby a simple rate law derived from classical models.It was found that there is an additional step in the dissolution process for aspirin, occurringbetween the initial rapid stage and the fixed rate law stage, during which the crystal attainsa solution annealed shape. The fixed rate law stage was again found to be described by asimple rate law.iiAbstractThe results obtained are in agreement with an earlier dissolution study of NaCl crystals,thus it appears that the classical rate laws can be used to describe the dissolution of a varietyof complex molecular and ionic crystals.iiiLay SummaryDissolution of molecular solids is an important process in many physical systems and situ-ations. However, a lot of work is still needed for a detailed understanding of the molecularmechanism of dissolution, because of the complex nature of many molecular solids, and thelarge time scales required for simulation studies. My study was focused on uncovering themolecular mechanism of dissolution, and employed molecular dynamics simulation to inves-tigate the complete dissolution of molecular solids. Urea and aspirin molecules are used asmodels. I found that these crystals dissolve in a simple three step process, and under sinkconditions, detachment of the molecules from crystal is the rate limiting step. I also foundthat the dissolution profile can be described by simple rate laws to a very good extent. Thus,it appears that these laws can be used to describe the dissolution of a variety of complexmolecular and ionic crystals.ivPrefaceThe research presented in the thesis are based on work done by me in Dr. G. N. Patey’sgroup. The chapters 3 and 4 are part of articles that will be published by A. Anand and G.N. Patey.The projects were designed by me with Dr. G. N. Patey, and are extensions of a previousproject, done by G. Lanaro, doctoral student in Dr. G. N. Patey’s group, which investigatedthe dissolution of ionic nanocrystals of NaCl.In all the projects, I performed the simulations, developed programs for data analysis,and formulated the hypothesis with suggestions, and guidance from Dr. G. N. Patey. Themanuscripts for the articles have been written by me with revisions and additions by Dr. G.N. Patey.vTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiiList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixList of Symbols and Abbreviations . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiv1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Foundation of dissolution research . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Recent trends in dissolution research . . . . . . . . . . . . . . . . . . . . . . . 32 Models and Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52.2 Water models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52.3 Urea model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62.4 Aspirin model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82.5 Algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8viTable of Contents2.5.1 Molecular dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92.5.2 Periodic boundary conditions . . . . . . . . . . . . . . . . . . . . . . . . 112.5.3 Temperature and pressure control . . . . . . . . . . . . . . . . . . . . . 132.5.4 Constraint algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . 163 Dissolution of Urea Nanocrystals . . . . . . . . . . . . . . . . . . . . . . . . . . 183.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 183.2 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 193.3 Models and methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 213.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 253.4.1 Dissolution profiles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 253.4.2 Temperature dependence and activation energies . . . . . . . . . . . . . 333.5 Summary and conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 374 Dissolution of Aspirin Nanocrystals . . . . . . . . . . . . . . . . . . . . . . . . 394.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 394.2 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 404.3 Models and methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 414.4 Results and discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 464.4.1 Stages of dissolution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 464.4.2 Temperature dependence and activation energy . . . . . . . . . . . . . 534.4.3 Model dependence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 544.5 Summary and conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 555 Summary and Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 575.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 575.2 Future directions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60viiList of Tables2.1 List of geometric, coulumbic and LJ parameters of water models . . . . . . . . 62.2 List of geometric and coulumbic parameters of urea molecules . . . . . . . . . . 73.1 The GAFF LJ parameters and partial charges for urea23 . . . . . . . . . . . . . 223.2 A summary of the numbers of urea and water molecules used, the shapes ofthe initial crystals, and the temperatures at which simulations were performed. 233.3 Goodness of fit parameters R2 obtained for different rate laws. Values of thebest fit are shown in bold. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 313.4 Rate constants (ns−1) for urea dissolution in different water models at 300,320, and 340 K. The rate constants for cubic and tablet-shaped crystals wereobtained from fits to the cube root and square root laws, respectively. . . . . . 323.5 Activation energies (kJ mol−1) for cubic crystal dissolution and urea diffusion. 354.1 A summary of the number of aspirin and water molecules used, shapes of theinitial crystals, and the temperatures at which simulations were performed . . . 444.2 Goodness of fit parameters R2 obtained for different rate laws . . . . . . . . . . 524.3 Diffusion coefficients of aspirin in solution. . . . . . . . . . . . . . . . . . . . . . 52viiiList of Figures2.1 Illustration of water models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62.2 Illustration of urea molcule and the unit cell . . . . . . . . . . . . . . . . . . . . 72.3 Illustration of aspirin molecule and the unit cell . . . . . . . . . . . . . . . . . . 82.4 A schematic of the MD algoritm . . . . . . . . . . . . . . . . . . . . . . . . . . 92.5 Periodic boundary conditions in two dimensions . . . . . . . . . . . . . . . . . . 113.1 Dissolution profiles obtained with the SPC water model for cubic crystals ofdifferent size (1024 and 432 molecules), plotted as the number of molecules inthe crystal vs time. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 223.2 Dissolution profiles for cubic crystals obtained with different water models.N(t) is the number of molecules in the crystal at time t. Note that results fromtwo simulation runs are shown for the SPC/E model. . . . . . . . . . . . . . . . 253.3 Simulations at different concentrations using the SPC water model. Curves Aand B are for urea mole fractions of 0.03428 and 0.02524 (calculated assumingcomplete dissolution of the entire crystal). Note that there is no significantconcentration effect. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 263.4 Configurational snapshots corresponding to different points in the dissolutionprofile of a cubic crystal in SPC water. Each red sphere indicates a urea molecule. 273.5 Configurational snapshots of the last stages in the dissolution profile of thecubic crystal in SPC water. Each red sphere indicates a urea molecule. . . . . . 28ixList of Figures3.6 Quantities obtained as described in the text for a cubic crystal in SPC water.All quantities are plotted as functions of the number of urea neighbors in thefirst coordination shell. The results were obtained over the time interval 1 - 3ns of a dissolution run. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 283.7 Configurational snapshots corresponding to different point in the dissolutionprofile of a tablet-shaped crystal in SPC water. Each red sphere indicates aurea molecule. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 293.8 Comparison of the dissolution profiles of cubic and tablet-shaped crystals inSPC water. fcry(t) is the fraction of molecules remaining in the crystal. . . . . 303.9 The radial density profile (molecules nm−3) about the center of a cubic crystalin SPC water, measured at 0.75 ns. The dashed horizontal line in the insetindicates the saturation density. . . . . . . . . . . . . . . . . . . . . . . . . . . . 323.10 Fits to the cube root law for a cubic crystal in SPC water at 300, 320, and 340K. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 333.11 Fit to the Arrhenius equation of the rate constants obtained for a cubic crytalin SPC water at 300, 320, and 340 K. . . . . . . . . . . . . . . . . . . . . . . . 343.12 Diffusion coefficients of urea (top panel) and water (bottom panel) in solution.Experimental trendlines for urea60 and water61 are included. . . . . . . . . . . 354.1 Schmatic representation of the cross section of a cubic crystal along the (010)crstal plane. The (100) and (001) planes lie perpendicular to the plane of thepaper. The oxygen atoms are represented in red, hydrogen atoms in gray, andcarbon atoms in black, respectively. . . . . . . . . . . . . . . . . . . . . . . . . . 424.2 a) Side views (curved surface), and b) Top views (flat surface) of the cylindricalcrystals. The crystals with (001), (010) and (100) faces as the flat surface areshown in i), ii) and iii), respectively. The oxygen atoms are represented in red,hydrogen atoms in gray, and carbon atoms in black, respectively. . . . . . . . . 43xList of Figures4.3 Dissolution profiles for the cubic crystal with different radii of the sphere usedfor removal of molecules in the solution. The profiles shown in the inset are forthe same crystals with the sphere radii reversed at 240 ns in the dissolution runs. 464.4 The dissolution profile for a cubic crystal displayed as the number of moleculesvs time. The profiles for a cubic crystal and a cylindrical crystal with the (100)face as its flat surface are shown in the inset. . . . . . . . . . . . . . . . . . . . 474.5 Snapshots of the (100) face corresponding to different points in the dissolutionprofile of a cubic crystal. Each red sphere indicates an aspirin molecule. . . . . 484.6 Dissolution profiles for three cylindrical crystals with different flat surfaces. . . 484.7 Snapshots corresponding to different points in the dissolution profile of a cylin-drical crystal with the (100) face as its flat surface. Each red sphere indicatesan aspirin molecule. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 494.8 Snapshots of the last stages in the dissolution profile of a cubic crystal. Eachred sphere indicates an aspirin molecule. . . . . . . . . . . . . . . . . . . . . . . 524.9 Fits to the square root law for a cubic crystal at temperatures of 300, 320, and340 K. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 534.10 Fits to the Arrhenius equation of the rate constants obtained for a cubic crystalat temperatures of 300, 320, and 340 K. . . . . . . . . . . . . . . . . . . . . . . 544.11 Fits to the square root law for a cubic crystal in TIP4P/2005 water at temper-atures of 300, 320, and 340 K. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55xiList of Symbols and AbbreviationskB Boltzmann Constantd Bond Lengthθ Bond Angleφ Dihedral Angleq Point ChargeK Kinetic EnergyV Potential Energyf ForceD Diffusion CoefficientEa Activation EnergyuC Coulombic PotentialuLJ Lennard-Jones Potentialσ Lennard-Jones Length Parameter Lennard-Jones Energy ParameterAPI Active Pharmaceutical IngredientsGROMACS GROningen Machine for Chemical SimulationsLINCS Linear Constraint SolverLJ Lennard-JonesMD Molecular DynamicsPME Particle Mesh EwaldPBC Periodic Boundary ConditionsxiiList of Symbols and AbbreviationsSPC Single Point Charge ModelSPC/E Extended Single Point Charge ModelTIP3P Three-site Transferable Intermolecular Potential ModelTIP4P/2005 Four-site Transferable Intermolecular Potential Model of 2005GAFF Generalized AMBER Force FieldAMBER Assisted Model Building with Energy RefinementCHARMM Chemistry at HARvard Macromolecular MechanicsMMFF Merlock Molecular Force FieldxiiiAcknowledgementsI would like to thank my supervisor Dr. Gren Patey, for all his help, encouragement andenthusiastic discussions throughout these last two years. I would also like to thank the labmembers, and all the friends I have met here, for all the fun during work and life in general.And finally, I would like to thank my family for all the support.xivChapter 1Introduction1.1 Foundation of dissolution researchThe study of dissolution processes has been developing since the late 1800’s, and the firstdissolution study was reported in the literature by Arthur A. Noyes and Willis R. Whitney in1897 in an article1 entitled “The rate of solution of solid substances in their own solutions”.This laid the foundation for dissolution research for the next century, where different disso-lution models based on various assumptions were suggested as physical explanations for thedissolution process.Noyes and Whitney suggested that the dissolution rate can be predicted with high accuracyby considering the dissolution process as a process of diffusion, where the molecules diffuseinto the solution from a thin diffusion layer formed around the solid. The mathematicalexpression of the law that they proposed is the famous Noyes-Whitney equationdcdt= K(cs − ct) , (1.1)where dcdt is the dissolution rate, K is the rate constant, cs is the solubility of the substance,and ct is the concentration of the dissolved substance at time t.They performed an experiment measuring the concentration of benzoic acid and leadchloride in water and calculated the K value at various time points. They found the valuescalculated to be very similar, confirming their hypothesis of a “diffusion controlled dissolu-tion”.In 1900 Erich Brunner and Stanislaus von Tolloczko published an article2 which considered11.1. Foundation of dissolution researchthe dependence of dissolution on the exposed surface area. They proposed an expression whichwas derived from Equation (1.1) by letting, K = K1Sdcdt= K1S(cs − ct) , (1.2)where S is the exposed surface area of the solid. The models presented above did not considerany time dependent changes in the surface area of the solid during dissolution. Hixson andCrowell3 addressed the fact that the surface of the dissolving solid changes with time andobtained an expression for dissolution rate,dMdt= −K ′St(cs − ct) , (1.3)where dM is the amount of the substance which dissolved in the time interval dt, K′is aconstant, St is the available surface area, cs is the solubility in the bulk fluid, and ct is theconcentration of the dissolved substance at time t.For cases when the change in concentration of the dissolved substance in the bulk fluid isnegligible (perfect sink conditions), “cs − ct” can be considered as a constant. Thus k′ and“cs − ct” can be combined to a new constant k′′ such that,dMdt= −K ′′St . (1.4)For a spherical particle with radius Rt at time t the surface area, St can be written asSt = 4piR2t , the volume of the particle vt as, vt =43piR3t , and since “mass ∝ V olume”, onecan express surface area as a function of mass in the following way: St = constant ×M23t .Thus the Equation (1.4) can be written asdMdt= −K ′′′M23t , (1.5)with K′′′being another constant and Mt being the remaining mass of the solid at time t. The21.2. Recent trends in dissolution researchintegrated form of the Equation (1.5) leads to the well known “Hixson-Crowell equation” orthe “cube root law”M130 −M13t = kt . (1.6)where M0 and Mt denote the mass of the solid at time t = 0 and at time t, respectively, andk is a constant.There are different explanations that have been proposed for the rate law expressions listedabove. The first which assumes that diffusion through the surface layer is the limiting step ofthe dissolution process, was put forward by Brunner4 and Nernst.5 Another explanation camefrom Lanaro and Patey,6 where they show that under sink conditions, the rate determiningstep is the detachment of the molecule from the crystal surface. Various other models havealso been proposed as an alternative explanation to the dissolution process. The “interfacialbarrier model” which considered interfacial transport rather than the diffusion through thelayer as the limiting step, first proposed by Wilderman7 and later considered by Zdanovskii,8has not been studied extensively. Another model proposed by Dankwerts9 considered thatconstantly renewed packets of solvent reach the surface and absorb solute particles, deliveringthem to the solution. Combinations of these models have also been studied.1.2 Recent trends in dissolution researchAll the advances in in-vitro dissolution studies were achieved before 1950 by physical chemistslaying the basic set of principles for the dissolution process. The second half of the 20th centurysaw a considerable amount of interest in dissolution studies again due to a couple of factors.First the relationship between dissolution rates and bioavailability10,11 started to develop inthe late 1950’s and by the early 21st century in vitro dissolution testing emerged as a veryimportant tool for the development and approval of safe and generic drug products.12,13 Thissparked considerable interest in the dissolution process of various crystalline solids.14The second factor was the large increase in computational power, as a result of whichmolecular dynamics simulation studies of the dissolution process started to gain popularity in31.2. Recent trends in dissolution researchthe scientific community. These studies were conducted to uncover the steps of dissolution atthe microscopic level, as the classical models discussed before presented very little informationabout the dissolution of nanocrystals. Piana and Gale15 studied the dissolution and growthof the [001] surface of urea in contact with water. They found that the single surface dissolvesquickly, forming a super saturated solution, which on cooling leads to rapid growth of thesurface. Gao and Olsen16 studied dissolution of the drug acetaminophen in water, and revealedthe importance of corners and edges in the initial stages of the dissolution process.Elts et. al.17 studied the dissolution of aspirin using both Molecular dynamics simulationand Kinetic Monte Carlo techniques, and investigated the possibility of predicting dissolutionrates for Active Pharmaceutical Ingredients (API) based on their molecular structures. How-ever, because of the complex structures and low solubilities of these API’s, the computationaleffort required for the study of the whole dissolution process of bigger molecules is not prac-tical. Lanaro and Patey6 investigated the dissolution of an ionic crystal, NaCl, in water, andanalyzed the dissolution rate using some of the classical dissolution models discussed before.They found that these models provide a very good description of the intermediate stage ofthe dissolution profile. Using this technique they could predict the dissolution rate of NaClin water by analyzing only a fraction of the full dissolution run.Thus it is of interest to test if the same technique could be used for molecular solids, asthis would help reduce the time scale required for the prediction of rates for APIs.4Chapter 2Models and Methods2.1 OverviewSimulation requires a model of the system that mimics the experimental properties of thesample for the conditions under investigation. Simulation also requires various methods tocontrol the temperature, pressure, and the molecular geometry (chemical bonds and bondangles) of the system during evolution. In this chapter we discuss the models that are usedfor the dissolution studies of two different molecular solids (aspirin and urea) in water, andexplain the various algorithm used during the simulation.2.2 Water modelsThe different water models used in the dissolution studies are: Simple Point Charge (SPC) wa-ter model,18 Extended Simple Point Charge (SPC/E) water model,19 three point TransferableIntermolecular Potential (TIP3P)20 water model, and the four point Transferable Intermolecu-lar Potential (TIP4P/2005)21 water model. The first three models (SPC, SPC/E, and TIP3P)have three interaction sites (one on each atom of the molecule), while the TIP4P/2005 watermodel has four interaction sites, with the fourth virtual site representing the lone pair placedat the bisector of the H−O−H angle. The molecular geometries of the different water modelare displayed in Figure 2.1 and the geometric parameters, Lennard-Jonnes (LJ) parametersand the charges on the interaction sites are summarized in Table 2.1. The water models usedare rigid and have slightly different Coulombic and LJ parameters. Different models are usedto study the influence of the solvent model on the dissolution profile of molecular solids.52.3. Urea modelFigure 2.1: Geometrical representation of the three and four site water models. The oxygenatoms are represented in red, hydrogen atoms in gray, and the virtual sites in pink, respec-tively.SPC SPC/E TIP3P TIP4P/2005dO−H (nm) 0.1000 0.1000 0.0957 0.0957dO−M (nm) n/a n/a n/a 0.0155θH−O−H (degree) 109.5 109.5 104.5 104.5qO (e) -0.8200 -0.8476 -0.8340 0.0qH (e) 0.4100 0.4238 0.4170 0.5564qM (e) n/a n/a n/a -1.1128σO (nm) 0.3016 0.3016 0.3151 0.3158O (kJ mol−1) 0.6500 0.6500 0.6364 0.7749Table 2.1: Molecular geometry, LJ parameters, and partial charges on interaction sites fordifferent water models. The oxygen atom, hydrogen atoms, and virtual sites are denoted byO, H, and M respectively. The distance between sites is denoted by d, the bond angle by θ,the charge by q, and σ and  are the LJ parameters.2.3 Urea modelThere are a number of atomistic models for urea that have been put forward for simulationstudies.22 However, the model that best reproduces the observed molecular properties of ureacrystal and solution is the improved generalized AMBER force field (GAFF)23,24 model. Themolecular geometry of urea is displayed in Figure 2.2 and the improved GAFF parametersare summarized in Table 2.2. The unit cell representation was taken from X-ray diffraction62.3. Urea modelstudies, which have been conducted by many authors,25–27 and is displayed in Figure 2.2.Figure 2.2: Geometrical representation of a single urea molecule and the unit cell of theurea crystal. The space group is P 4¯21m and the unit cell dimensions are a = b = 0.5565and c = 4.684 nm. The oxygen atom is represented in red, hydrogen atoms in gray, nitrogenatoms in blue, and the carbon atom in black.dC−O (nm) dC−N (nm) dN−H (nm)0.1250 0.1383 0.1010KdC−O KdC−N KdN−H2744.70 1774.02 1815.86θN−C−O (degree) θC−N−H (degree) θH−N−H (degree) θN−C−N (degree)120.9 120.0 120.0 118.6KθN−C−O KθC−N−H KθH−N−H KθN−C−N334.72 125.52 146.44 292.88φH−N−C−O (degree) φH−N−C−O (degree) φN−N−C−O (degree) φC−H−N−H (degree)180.0 0.0 180.0 180.0KφH−N−C−O KφH−N−C−O KφN−N−C−O KφC−H−N−H10.46 8.37 43.93 4.60qO (e) qC (e) qN (e) qH (e)-0.660 0.884 -0.888 0.388Table 2.2: Improved GAFF parameters for urea.23 The oxygen atom, carbon atom, hydrogenatoms, and nitrogen atoms are denoted by O, C, H, and N, respectively. The bond lengthis denoted by d, bond angle by θ, dihedral angle by φ, and the charges on the sites by q.The force constants, Kd are given in kJ mol−1 A˚−2, Kθ in kJ mol−1 radian−2, and Kφ in kJmol−1.72.4. Aspirin model2.4 Aspirin modelThe unit cell representation of the aspirin crystal has been taken from a neutron diffractionstudy.28 The unit cell was relaxed for the force field, and the modified lattice parameters29obtained from optimization of the bulk crystal structure were in good agreement with theexperimental data. The molecular geometry of the aspirin molecule and the unit cell aredisplayed in Figure 2.3.Figure 2.3: Geometrical representation of the aspirin molecule and the unit cell of aspirincrystal. The space group is P21/c. The modified lattice parameters are, a = 1.140, b = 0.660,c = 1.150 nm, and β = 91.9 degrees. The oxygen atoms are represented in red, hydrogenatoms in gray, and the carbon atoms in black.2.5 AlgorithmsIn this section some of the methods used by GROMACS during the simulation of varioussystems are elucidated. First, the Molecular Dynamics(MD) algorithm is discussed, followedby the periodic boundary conditions, temperature and pressure coupling methods, and finally82.5. Algorithmsthe constraint algorithm is discussed in detail.2.5.1 Molecular dynamicsFigure 2.4: The basic Molecular Dynamics (MD) algorithmMolecular dynamics is a computational method used to study the evolution of a N-bodysystem. A schematic of the basic MD algorithm is displayed in Figure 2.4. The initialcoordinates and the velocities (can also be generated by GROMACS) are specified as input,along with a description of the interatomic interactions and the system conditions under study.The initial states for simulation were generated using the models discussed in the Sections2.2, 2.3, and 2.4. Initial configurations were obtained by placing the crystal at the center ofa cubic simulation cell and the filling the cell with water molecules. The system is allowed toevolve in time and the trajectories are obtained by numerically solving Newton’s equations ofmotion,fi = mir¨i , (2.1)where mi is the mass of the corresponding interaction site, r¨i is the acceleration, and fi is thetotal force experienced by site i. The trajectory, ri(t), is obtained by numerically integrating92.5. AlgorithmsNewton’s equations of motion, and can be used to calculate various dynamical and equilibriumproperties of the system under investigation. The total force fi on any site is computed byconsidering both bonded and non-bonded interactions plus any restraint or external force.The method used to integrate the equations of motion is the leap-frog algorithm,30 whichupdates the velocities and positions of the particles as follows:ri(t+ ∆t) = ri(t) + ∆tvi(t+∆t2) , (2.2)vi(t+∆t2) = vi(t− ∆t2) + ∆tai(t) , (2.3)ai(t) =fi(t)mi, (2.4)fi(t) = −∇ui(ri) . (2.5)It can be noted that the positions and velocities are calculated at different time steps. A verysmall time step, ∆t (1 femtosecond) has been chosen for all the simulations performed in thisthesis.The bonded interaction parameters (force constants) are included in the force field modeland the bonds are constrained during the evolution of the system using the method describedin section 2.5.4. The non-bonded interaction are modeled using the LJ and coulombic poten-tials, and are computed by considering a list of non-bonded atoms within a specified cutoffradius. The LJ potential consists of two terms, the r−12 term which is a repulsive termdescribing the short range repulsion due to overlapping electron orbitals, and the r−6 termwhich is an attractive term describing attraction at long rangeuLJ(rij) = 4ij[(σijrij)12−(σijrij)6], (2.6)where σij and ij are the LJ parameters for the sites under consideration, and rij is the inter-atomic distance. The Lorentz-Berthelot combining rules31 are most often used to calculate102.5. Algorithmsthe LJ parameters for a pair of atoms from individual atom LJ parameters,σij =12(σi + σj) , (2.7)ij = (ij)12 , (2.8)where σi, σj , i, and j are the individual particle LJ parameters. The coulombic potentialdepends on the partial charges on the atoms, qi and qj , and the interatomic separation rij ,uC(rij) =14pi0qiqjrij, (2.9)where 0 is the permittivity of free space.2.5.2 Periodic boundary conditionsFigure 2.5: Periodic boundary conditions in two dimensions. The dark blue box is the centralcell which is repeated in all directions.There are various conditions that are employed to calculate bulk properties by using a small(finite) part of the whole (infinite) system. These conditions known as the periodic boundary112.5. Algorithmsconditions (PBC), minimize the edge effects of the finite system under investigation. ThePBC algorithm surrounds the finite box by translated copies of itself, as shown in Figure 2.5,and as the system evolves particles in all the boxes move in an identical manner. Thus, if aparticle leaves the box at one end, an identical particle enters the box at the other end.As a result of PBC, there is an enormous increase in the number of interacting pairs, asthe particle not only interacts with other particles in the box but with their images. Thisproblem is tackled by using a minimum image convention combined with a spherical cutoff,which only allows the nearest neighbors of particle images to interact.The short-range interactions in the minimum image convention are calculated by consid-ering the interaction of an atom with the nearest atom or image in the periodic system. Ifan atom leaves the box, the interaction calculation is done by considering the incoming im-age. The cutoff for short-range interactions is restricted to half the box length, which can besignificantly bigger than the cutoff for the potential used (of the order of ∼1 nm).The long-range electrostatic forces have ranges greater than half the box length and arethus harder to treat. The total electrostatic energy of an atom in the central box can bewritten asui =14pi0∑nN∑j′ qiqjrij,n, (2.10)where n = (nx, ny, nz) is the box index vector, rij,n is the actual distance between the chargesand the prime indicates that i 6= j when n = {0}. The sum in Equation (2.10) is conditionallyconvergent which is troublesome and a number of methods have been proposed to tackle thisproblem.The Ewald summation32 is a technique that is used to calculate the long-range interactions.It basically converts the single long-ranged term into two short-ranged terms (one in real spaceand the other in reciprocal space) and a constant term. The expression for the electrostaticenergy after using Ewald summation can be written asui = uSi + uLi + uselfi , (2.11)122.5. AlgorithmsuSi =14pi012∑nN∑j′ qiqjrij,nerfc(rij,n212σ), (2.12)uLi =12V 0∑k 6=0N∑jqiqjk2expιk.rij exp−σ2k2/2 , (2.13)uselfi = −14pi01(2pi)12σq2i , (2.14)where erf(z) = 2pi1/2∫ z0 exp−t2 dt, erfc(z) = 1 − erf(z), V is the volume of the cell, k = |k|, kis the reciprocal lattice vector, and σ is the standard deviation of the Gaussian distributionwhich is used to modulate the convergence rate. The first term in Equation (2.11), uSi issimilar to the total interaction energy ui, and can be truncated at a certain cutoff due toinclusion of the erfc(r)/r term, which decays very rapidly with r. The second term uLi is thesum of the long-range potential in reciprocal space, and uselfi is the self interaction term, tonegate the extra term added in the second term uLi , to include the potential generated by thesite i itself.The classical Ewald method scales as O(N2), while the variant particle mesh Ewald33(PME) method is an advanced method using the fast Fourier transform34 technique to com-pute the reciprocal term, which scales as O(N log(N)), where N is the number of charges.Thus the PME method is most commonly used as it is more efficient.2.5.3 Temperature and pressure controlThe molecular dynamics algorithm discussed in Section 2.5.1 gives rise to the NV E (constantnumber of particles, volume, and energy) ensemble. However, most of the quantities wewish to calculate are at fixed temperature and pressure, which cannot be well controlled byusing the above algorithm. So, we now discuss various methods that are used to control thetemperature and pressure of a system during a simulation.Nose´-Hoover thermostatThe instantaneous temperature of a system of N particles is related to the kinetic energy132.5. Algorithmsof the system K byT =2KkBNf, (2.15)where Nf is the total number of degrees of freedom. The Nose´-Hoover thermostat35,36 in-troduces terms representing a thermal reservoir into the system Hamiltonian and a dynamicparameter ξ, whose physical meaning is that of a friction which accelerates or decelerates theparticles until the temperature is equal to the desired value. The final Hamiltonian of thesystem after inclusion of the parameter isH = K − V + p2ξ2Q+NfkBT0ξ , (2.16)where pξ is the conjugate momentum of the friction parameter, and K and V are the kineticand potential energies of the system of N particles, respectively. The last two terms inEquation (2.16) can be interpreted as the kinetic and potential energy associated with athermal reservoir. The equations of motion are modified tor¨i =Fimi− pξQr˙i , (2.17)dpξdt=N∑imiv2i2− Nf + 12kBT0 , (2.18)where T0 is the target temperature and Q is a mass parameter which along with the targettemperature determines the strength of the coupling. The mass parameter Q is related to thetarget temperature T0 byQ =τ2TT04pi2, (2.19)where τT is the period of oscillation of kinetic energy between the system and the reservoir,and is a better parameter to control the coupling strength, because it is independent of thesystem size and target temperature.Parinello-Rahman barostatThe Parinello-Rahman barostat37 like the Nose´-Hoover thermostat modifies the equations142.5. Algorithmsof motion to allow the box vectors to evolve in time to produce the desired pressure. The boxvectors a, b and c are time dependent and are represented by a matrix h = {a, b, c}. Thevolume of the box V is then defined as: V = det(h) = a • (b× c), where the dot representsa scalar product. The Lagrangian of the system subjected to an external pressure p can bewritten asL = K − V + 12MTr(h˙′h˙)− pV , (2.20)where the prime denotes the transpose, and K and V are the kinetic and potential energyof the system of N particles, respectively. The last two terms in Equation (2.20) can beinterpreted as the kinetic and potential energies associated with a piston applying an externalpressure. The coordinates of the particles are now represented assi = h−1ri , (2.21)and the equation of motion is modified tos¨i = −∑j 6=idudrij1mirij(si− sj)−G−1G˙s˙i , (2.22)h¨ =1M(pi − p)σ , (2.23)where the terms pi, the microscopic stress tensor, σ, a matrix containing the direction of thereciprocal vectors and G, known as the metric stress tensor are defined as follows:V pi =∑imiviv′i −∑i∑j>idudrij1rijrijr′ji , (2.24)σ = {b× c, c× a, a× b} , (2.25)G = h′h . (2.26)152.5. Algorithms2.5.4 Constraint algorithmIn a molecular dynamics simulation, the sites move around under the influence of intermolec-ular forces affecting the geometry (chemical bonds and bond angles) of some molecules. Thuswe need to apply certain constraints on every molecule to conserve the molecular geometryduring the evolution. In this section we discuss the linear constraint solver (LINCS)38 al-gorithm which is a method that resets the bonds to their correct lengths in two steps, anunconstrained evolution and a projection of the new bonds on the old ones.The system is constrained using a set of k time independent equationsgi(r) = |ri1 − ri2 | − d = 0, i = 1, 2, ....k , (2.27)r = {r1, r2.....rN} , (2.28)where d is the length of the bond between atoms i1 and i2. These holonomic constraints(dependent only on position r and time t) are included in the Lagrangian asL = K − V + λ • g , (2.29)where λ = {λ1, ...λk} are the k Lagrange multipliers associated with the constraints g ={g1, ...gk}, and K and V are the kinetic and potential energies of the system. The resultingequations of motion can be expressed asr¨ = (I− TB)M−1f −TB˙r˙ , (2.30)T = M−1B′(BM−1B′) , (2.31)where the prime and dot represent the transpose and derivative with respect to time, respec-tively, M = diag(m1,m1,m1, ....mN ,mN ,mN ) is a 3N dimensional diagonal matrix contain-ing the masses of the particles, and B is a K × 3N matrix, containing the directions of the162.5. AlgorithmsconstraintsBhi =∂gh∂ri. (2.32)I− TB is a projection matrix, M−1f is the acceleration vector, and the last term in Equation(2.30) represents centripetal forces caused by rotation of bonds. The LINCS algorithm can beeasily implemented by first using the leap-frog algorithm to find the unconstrained positionsand velocities (r∗ and v∗), and then projecting them to conserve the constraints as follows:r(t+ ∆t) = (I− TB)r∗(t+ ∆t) + Td , (2.33)v(t+∆t2) = (I− TB)v∗(t+ ∆t2) +1∆tT(Br− d) , (2.34)where the terms Td and 1∆tT(Br− d) are position and velocity correction terms added tothe equations to prevent accumulation of numerical error in computer implementation,38 andd is a vector containing the bond lengths.17Chapter 3Dissolution of Urea Nanocrystals3.1 OverviewMolecular dynamics simulations are used to determine the mechanism of urea crystal dis-solution in water under sink conditions. Crystals of cubic and tablet shape are considered,and results are reported for four commonly used water models. The dissolution rates for thedifferent water models differ considerably, but the overall dissolution mechanism remains thesame. Urea dissolution occurs in three stages: a relatively fast initial stage, where moleculesleave high energy sites such as corners and edges, a slower intermediate stage during whichthe bulk of the crystal dissolves, and a final stage where the small crystal remnant loosesits structure and completely disintegrates. We show that during the long intermediate stagethe dissolution process is well described by a simple classical model which assumes that thedissolution rate is proportional to the active surface area. The active surface area is a regionof the crystal surface from which molecules leave uniformly after the initial stage of dissolu-tion. For initially cubic and tablet-shaped crystals the active surface areas are spherical andcylindrical, respectively, and the corresponding integrated rate laws give excellent fits to thesimulation results. By carrying out simulations at different temperatures we show that ureadissolution is an activated process, the rate constants obey the Arrhenius equation with anactivation energy of ∼ 32 kJ mol−1 for a cubic crystal. Our simulations give no indicationof a significant diffusion layer, and we conclude that the detachment of urea molecules fromthe crystal surface is the rate determining step in the dissolution process. We note that thedissolution mechanism, and the applicability of classical rate laws, that we report for urea areconsistent with earlier observations for the dissolution of NaCl crystals.6 This suggests that183.2. Introductionthe three-stage mechanism and classical rate laws might also apply to the dissolution of otherionic and molecular crystals.3.2 IntroductionMolecular solids are an important class of compounds of relevance in the pharmaceutical andspecial chemicals industries, among other applications.12,13,39–44 The dissolution of solids isa crucial process in many physical systems and situations.12,13 Given this, crystal dissolutionhas been a subject of interest for the past century,10,45 with current research strongly focussedon gaining a detailed understanding of dissolution at the molecular level.14In principle, computer simulations employing molecular dynamics (MD) or other methodsprovide an ideal approach to gaining a microscopic understanding of dissolution. However,in practice, until recently this method has not been widely exploited due partially to thecomplex nature of interesting molecules, and because the large amounts of simulation timerequired exceeded the power of most computational facilities. However, during the past fewyears simulations of dissolution have become feasible due to great increases in computingpower, and to the ever evolving development of accurate forcefields for complex molecules.We note that the last point is particularly important in the dissolution context, since to beuseful in dissolution studies a model must give a good representation of both solution andsolid properties.Despite the potential issues, there have been several recent investigations of dissolutionemploying simulation methods.6,15–17 We note in particular the work of Gao and Olsen16which revealed the importance of edges and corners as the sites from which molecules initiallydetach from a crystal. In another important study, Elts et al.17 have proposed a simulationapproach that combines MD and kinetic Monte Carlo methods, such as to overcome the longsimulation times required for the dissolution of some crystals. They applied their techniquesin a detailed study of aspirin dissolution in water, and reported a dissolution mechanism ingood accord with experimental observation.193.2. IntroductionIn a recent paper, Patey and Lanaro6 employed direct MD simulations to investigate thecomplete dissolution of differently shaped NaCl nanocrystals. For model NaCl, dissolutionis sufficiently fast that the dissolution process can be followed until an entire crystal hasdissolved. An interesting finding of this work was that NaCl dissolution essentially occurs inthree stages, a relatively fast initial stage where particles are removed from high energy sites(mainly edges and corners, depending on the initial crystal shape), a long slower intermediatestage where dissolution closely follows a fixed rate law, and a final stage where crystallineorder is lost as the dissolution process is completed. Additionally, it was shown that duringthe intermediate stage which dominates the dissolution process, the dissolution rate is welldescribed by classical models,1–5,46,47 which take account of the varying surface area of thecrystal as it dissolves. Ion detachment from the NaCl crystal is the rate determining step,and the particular form of the integrated rate law depends on the initial crystal shape.The purpose of the present chapter is to examine if, and to what extent, the rather simplethree-stage dissolution process found for NaCl holds for molecular solids. To do this, weexamine the complete dissolution of urea nanocrystals in water under so-called sink conditions,where the solution concentration at complete dissolution is much less than saturation. Wenote that there have been many simulation studies of model aqueous urea solutions,15,48–50and that Piana and Gale15 have examined both crystallization and dissolution of urea inaqueous systems. However, while their work is interesting, Piana and Gale only consideredsystems where a single urea surface (001) was in contact with water, and do not address thestages of dissolution, and the shape-dependent rate laws in which we are interested.We selected urea for this investigation for several reasons. Firstly, as noted above, a studyof this type requires a forcefield which gives a good description of both crystalline and solutionphases. As discussed in Section 3.3, this is not true of all common urea models, but there isat least one urea model which satisfies both conditions. Secondly, at ambient temperatures,the dissolution of urea occurs sufficiently fast that nanocrystals (∼ 1000 urea molecules) canbe completely dissolved on feasible simulation time scales, allowing us to follow the completedissolution process. Finally, urea is a molecule of importance in many physical systems.51–53203.3. Models and methodsIn this chapter, we examine the dissolution of both cubic and tablet-shaped urea nanocrys-tals at temperatures ranging from 300 to 340 K. Four different water models are considered,and, although the observed dissolution processes are qualitatively similar, some significantquantitative differences do occur. Our main finding is that the dissolution of urea crystalsfollows a pattern that is very similar to that observed for NaCl, with simple rates laws givinga good description of the important intermediate stage of dissolution. Our analysis of ratelaws and the activation energy indicates that detachment of urea molecules from the crystalis the rate determining step in the dissolution process.The remainder of the chapter is organized into three sections. The models and methodsare described in Section 3.3, a detailed description of our findings is given in Section 3.4, andour conclusions are summarized in Section Models and methodsIn the present simulations all site-site pair interactions, u(rij), consist of Lennard-Jones (LJ)plus Coulombic terms such thatu(rij) = 4ij[(σijrij)12−(σijrij)6]+14pi0qiqjrij, (3.1)where qi and qj are charges on sites i and j, rij is the distance between sites i and j, σij andij are the LJ length and energy parameters, and 0 is the permittivity of free space.During the past few decades a number of different atomistic models for urea have beenproposed.22 For example, the OPLS/GROMOS54,55 model for urea has been widely employedin simulations of aqueous solutions. However, we found that with this model urea nanocrystalswere not stable under ambient conditions.In order to investigate urea dissolution, we need a model that produces stable nanocrystalsof urea. The best model available for this purpose appears to be the the Generalised AMBERForce Field (GAFF),23 which has been shown to give stable urea nanocrystals, and a gooddescription of the properties of crystalline urea.24 The GAFF potential parameters are given213.3. Models and methodsin Table 3.1. The GAFF parameters are used to generate an AMBER topology file, which isthen converted to GROMACS topology for simulations. Several different water models havebeen used in studies involving urea in aqueous solution.15,22,50 Among others, these includethe SPC,18 SPC/E,19 TIP3P,20 and TIP4P56 water models. In the present chapter we carryout simulations with the SPC, SPC/E, TIP3P, and TIP4P/200521 models to investigate anydependence on the particular water model employed. We note that, of all the rigid watermodels available, TIP4P/2005 gives the best overall agreement with experimental results forpure water.57 For all models the usual Lorenz-Berthelot combining rules31 were used in thecalculation of the LJ interactions.σ(A˚) (kJ mol-1) q(e)C 1.9080 0.3598 0.884O 1.6612 0.8786 -0.660N 1.8240 0.7113 -0.888H 0.6000 0.0657 0.388Table 3.1: The GAFF LJ parameters and partial charges for urea23Figure 3.1: Dissolution profiles obtained with the SPC water model for cubic crystals ofdifferent size (1024 and 432 molecules), plotted as the number of molecules in the crystal vstime.223.3. Models and methodsSimulations were carried out for urea nanocrystals shaped as cubes and as tablets, in orderto examine the influence of shape on the dissolution process. The cubic crystals consistedof 1024 urea molecules, and were constructed by repeating the urea unit cell26 eight timesin each direction. In order to ensure that our conclusions on dissolution are not stronglydependent on the crystal size, a test simulation using SPC water was performed with a crystalof 432 molecules at 300 K. The results were qualitatively similar to those obtained with 1024molecules as can be seen in the Figure 3.1. Tablet-shaped crystals contained 1264 molecules,and were generated by repeating the unit cell 16, 16, and 4 times along the crystal axes, andthen removing all molecules that were outside a circle of 4 nm centered on the axis of thecylinder. Configurational snapshots of both cubic and tablet-shaped crystals are given below(see Figures 3.4 and 3.7).The experimental saturation concentrations48 of urea for our simulation conditions cor-respond to urea mole fractions in the range 0.26 − 0.43. All simulations were carried out aturea mole fractions that were less than 0.035, or at concentrations less than one eighth ofthe experimental saturation concentrations. We would expect these concentrations to meetthe so-called sink condition, which refers to solutions where the concentration is so low thatmolecular reattachments do not significantly influence the dissolution rate.58 We tested thatsink conditions were indeed met for our model system by carrying out simulations at differentconcentrations, as discussed in Section 3.4. The numbers of urea and water molecules includedin each simulation, together with the temperatures considered, are summarized in Table 3.2.Urea Water Models Temp(K)SPC SPC/E TIP3P TIP4P/2005Cube 1024 29873 29873 29873 30125 3001024 29873 29873 29873 30125 3201024 29873 29873 29873 30125 3401024 40576 n/a n/a n/a 300Tablet 1264 39672 n/a n/a n/a 300Table 3.2: A summary of the numbers of urea and water molecules used, the shapes of theinitial crystals, and the temperatures at which simulations were performed.All simulation were carried out under NPT conditions at a pressure of 1 bar, using the233.3. Models and methodsGROMACS59 molecular dynamics package, version 4.5.5. The temperature was controlledusing the Nose´-Hoover thermostat35,36 with a time constant of 2.0 ps, and the pressure wasmaintained by means of the Parrinello-Rahman barostat,37 with a compressibility of 4.46 ×10−5 bar−1 and a time constant of 2.0 ps. The time step used for all simulations was 1 fs.Periodic boundary conditions were applied, and the short-range interactions were truncatedat 1 nm. The electrostatic interactions were taken into account using the fourth order particlemesh Ewald method33 with a Fourier spacing of 0.12 nm. All bonds were constrained usingthe LINCS algorithm.38Initial configurations were obtained by placing the urea crystal at the center of the cubicsimulation cell, and the filling the cell with water molecules. Water molecules located too closeto the surface of the crystal (within 0.1 nm) were removed. The system was then relaxed fora brief period, keeping the urea molecules fixed, so as to remove any gaps occurring at thecrystal-water interface due to the removal of water molecules. The system was then evolvedin time until complete dissolution of the crystal. This process took times ranging from 2.10ns to 12.05 ns, depending on the temperature of the system.In order to identify urea molecules as solution or crystal molecules, an order parameterbased on counting the number of neighboring urea molecules was used. The number ofneighbors of a molecule in the crystal is substantially higher than that of a molecule insolution, and this difference was used to devise a simple order parameter for classification.For this purpose, the number of neighbors is defined as the number of molecules within asphere of radius 0.6 nm, centered at the center of mass of the molecule. The number of ureaneighbors of a molecule in the crystal lies in the range 12− 14, and that of a molecule in thesolution is in the range 0− 2. Therefore, we chose 7 neighbors as a cutoff value, any moleculewith fewer than 7 neighbors is considered to be part of the solution, and those with 7 ormore neighbors are considered to be part of the crystal. The simulation results were not verysensitive to the exact value of the cutoff, similar dissolution profiles were for cutoffs rangingfrom 6 to 8 neighbors.243.4. Discussion3.4 Discussion3.4.1 Dissolution profilesCubic crystalsDissolution profiles for an initially cubic urea crystal at 300 K are plotted in Figure 3.2.Results are included for all four water models, and we see some significant model dependence.The time required for complete dissolution varies from ∼ 7 ns in the fastest case (TIP3P) to∼ 11 ns in the slowest (SPC/E), with SPC and TIP4P/2005 both taking ∼ 8 ns. For all fourmodels the initial conditions were very similar, and consistent model differences occurred indifferent simulation runs, and at different temperatures. Note that profiles for two simulationswith SPC/E water with different initial conditions are included in Figure 3.2, and that onlysmall differences are observed. Therefore, the discrepancies apparent in Figure 3.2 reflect areal dependence on the water model employed. A further discussion of the model dependenceis given below, where we consider the temperature dependence and activation energies (Section3.4.2).Figure 3.2: Dissolution profiles for cubic crystals obtained with different water models. N(t)is the number of molecules in the crystal at time t. Note that results from two simulationruns are shown for the SPC/E model.253.4. DiscussionIn order to be certain that we are truly under sink conditions and that finite solutionconcentration is not influencing our results, we carried out two simulations (using SPC water)at mole fractions of 0.03428 and 0.02561 (calculated at complete dissolution), keeping allother conditions fixed. The dissolution profiles are shown in Figure 3.3, and we see that nosignificant differences are observed, confirming that we are indeed in the sink regime.Figure 3.3: Simulations at different concentrations using the SPC water model. Curves A andB are for urea mole fractions of 0.03428 and 0.02524 (calculated assuming complete dissolutionof the entire crystal). Note that there is no significant concentration effect.Despite the water model dependence of the overall dissolution rate, by closely examiningdissolution profiles, together with configurational snapshots, we determine that for all watermodels the urea dissolution process can be divided into three stages, that are similar to thosepreviously observed6 for the ionic crystal NaCl. An initial stage, where molecules are detachedfrom the edges and corners of the crystal, followed by an intermediate regime where the crystaldissolves according to a fixed rate law, after it has lost its corners and edges, and a final stagewhere the solid looses its crystalline structure and completes its dissolution. The following isa more detailed discussion of the three dissolution stages.During the initial stage of dissolution water molecules interact with the loosely boundmolecules generally located on edges and corners of the urea crystal, and these molecules are263.4. Discussionquickly dissolved. Following these rapid detachments, the active surface of the crystal furtheranneals such as to reduce the surface contact area for a crystal of given size. Configurationalsnapshots of a cubic crystal at various points on the dissolution curve are shown in Figure3.4. The snapshots shown were obtained with the SPC water model, but the particular watermodel used does not change the qualitative nature of the dissolution process.Figure 3.4: Configurational snapshots corresponding to different points in the dissolutionprofile of a cubic crystal in SPC water. Each red sphere indicates a urea molecule.From the dissolution profiles shown in Figure 3.2 (note the magnification of the initialstage), we see that the initial slopes (between 0 and ∼ 0.25 ns) are quite steep, after whichthe slopes decrease and remain roughly constant throughout the stage we refer to as the fixedrate law regime. During the rapid initial stage ∼ 150 urea molecules (∼ 14%) leave the crystaland enter the solution. It is evident from the snapshots in Figure 3.4 that these moleculesleave the edges and corners of the crystal. After the initial period more molecules detach fromthe crystal, and this continues until the crystal is nearly spherical in shape, as can be seenin Figure 3.4(d). Once the crystal attains a spherical shape, the molecules leave the surfaceessentially uniformly, thus entering the fixed rate law regime. This regime continues until thefinal stage of dissolution (approximately the last nanosecond of the profiles shown in Figure3.2), when the solid loses its crystalline order, becoming an amorphous urea cluster whichcontinues to disperse into solution until no crystal remnant remains. This is illustrated in theconfigurational snapshots shown in Figure DiscussionFigure 3.5: Configurational snapshots of the last stages in the dissolution profile of the cubiccrystal in SPC water. Each red sphere indicates a urea molecule.Figure 3.6: Quantities obtained as described in the text for a cubic crystal in SPC water. Allquantities are plotted as functions of the number of urea neighbors in the first coordinationshell. The results were obtained over the time interval 1 - 3 ns of a dissolution run.To gain further insight into the influence of local environment on urea detachment fromthe crystal, we calculated the probability of detachment during a given time interval of surfacemolecules with different numbers of urea molecules in the first coordination shell. The firstcoordination shell is defined by a sphere of radius 0.5 nm centered at the molecular center ofmass. The radius of the first coordination shell was selected based on the urea carbon-carbonradial distribution function. Surface molecules were identified as crystal molecules (definedas described in Section 3.3) having 1− 6 urea molecules in their first coordination shell. Thenumber of molecules that detached from the crystal within 0.1 ns was recorded, and averages283.4. Discussionwere taken over 40 equally spaced time slices between 1 and 3 ns of a dissolution run. Thecalculation was started at 1 ns to ensure that the results were not influenced by the initialstructure of the surface. The results obtained with SPC water are shown in Figure 3.6, andwe see that the probability of detachment decreases dramatically as the number of neighborsincreases, as one might expect. We note that the the particular time interval considered isnot an important factor, results obtained using longer time intervals (0.15 ns and 0.2 ns) werequalitatively similar to those shown in Figure 3.6.Figure 3.7: Configurational snapshots corresponding to different point in the dissolution profileof a tablet-shaped crystal in SPC water. Each red sphere indicates a urea molecule.Influence of crystal shapeIn order to investigate the possible influence of crystal shape, we also consider a tablet-shapedcrystal, as depicted in Figure 3.7. A dissolution profile for the tablet obtained with SPC wateris shown in Figure 3.8, and compared with the corresponding result for a cubic crystal. Sincethe numbers of molecules initially in the tablet-shaped and cubic crystals differ a little (Table3.2), in Figure 3.8 we plot the fraction of molecules in the crystal, fcry(t) = N(t)/N0, whereN(t) is the number of molecules in the crystal at time t and N0 is the initial number ofmolecules in the crystal. From Figure 3.8 we see that both tablet-shaped and cubic crystalshave similar dissolution profiles, with an initial steep slope followed by a nearly constantdissolution rate. The initial slope for the tablet-shaped crystal can be explained by thefact that the weakly attached molecules at the edges of the tablet leave very quickly. Also,293.4. Discussionmolecules do not leave from the flat surface of the tablet, but only from the curved surface,such that the cylindrical shape of the crystal remains intact throughout the dissolution, as canbe seen from the configurational snapshots in Figure 3.7. After the initial stage, the tabletdissolves at a faster rate than the cubic crystal. This can be attributed to the fact that thenumber of layers in the tablet is half the number in the cubic crystal. Thus, molecules onthe active surface (the sides) of the tablet interact more weakly with the bulk crystal thanmolecules on the active surface (spherical) of the initially cubic crystal.Figure 3.8: Comparison of the dissolution profiles of cubic and tablet-shaped crystals in SPCwater. fcry(t) is the fraction of molecules remaining in the crystal.Comparison with dissolution modelsAs noted above, dissolution processes have been studied for over a century, and several possiblerate laws have been proposed.10,14 In earlier work6 on the dissolution of NaCl crystals, weshowed that a relatively simple rate law does indeed give a good description of dissolution inthe fixed rate law regime, and it is interesting to check if this is also the case for urea. Undersink conditions we might expect the detachment process to be rate determining, such thatthe dissolution rate could be described by303.4. DiscussiondN(t)dt= −kSactive(t) , (3.2)where N(t) is the number of molecules in the crystal at time t, Sactive(t) is the active surfacearea from which detachments occur, and k is a rate constant. This law, which takes accountof the changing surface area, was first suggested by Brunner and Tolloczko.2Integrated rate laws can be obtained by identifying different forms for Sactive(t) and in-tegrating eq 3.2. Noting that the active surface areas in the fixed rate law regime are ap-proximately spherical and cylindrical for the cube and tablet, respectively, functional formsfor Sactive(t) can be obtained by considering surface area to volume relationships. This givesSactive(t) ∝ N2/3(t) for the cube and Sactive(t) ∝ N1/2(t) for the tablet, and the correspondingintegrated rate laws3√N(t) = 3√N0 − kt, (cube), (3.3)√N(t) =√N0 − kt, (tablet), (3.4)where N0 is the number of molecules in the crystal at t = 0, and k will obviously be differentfor each rate law.Cube TabletSPC SPC/E TIP3P TIP4P/2005linear 0.98839 0.98244 0.98227 0.99435 0.98561cube root 0.99068 0.99339 0.99421 0.98790 0.99306square root 0.99312 0.99326 0.99381 0.99294 0.99839Table 3.3: Goodness of fit parameters R2 obtained for different rate laws. Values of the bestfit are shown in bold.We examined the suggested rate laws by fitting to the dissolution profiles, excluding theinitial and final parts of the profile which lie outside the fixed rate law regime. Goodness offit parameters R2 are given in Table 3.3. Results for a linear rate law are also included. FromTable 3.3, we see that the square root law gives the best fit for the tablet, as expected. Forthe cubic crystal, the expected cube root law gives the best fit for the SPC/E and TIP3P313.4. Discussionmodels, but not for SPC and TIP4P/2005, for which better fits are given by the square andlinear laws, respectively. However, we note that in general all three rates laws give good fitsto the dissolution profiles, and the differences in the R2 values are too small to clearly identifythe “best” rate law. One would need profiles on longer time scales (greater than the timeneeded for complete dissolution for the urea nanocrystals) to make clear distinctions. Therate constants obtained using the cube root law for cubic crystals and the square root law forthe tablet are given in Table 3.4. These rate constants clearly show the model dependencenoted above.CubeTemperature (K) 300 320 340SPC 0.57 ± 0.15 1.22 ± 0.12 2.71 ± 0.35SPC/E 0.40 ± 0.08 0.88 ± 0.22 1.84 ± 0.36TIP3P 0.65 ± 0.11 1.24 ± 0.20 2.63 ± 0.35TIP4P/2005 0.58 ± 0.16 1.52 ± 0.21 2.64 ± 0.23TabletSPC 6.26 ± 0.03 n/a n/aTable 3.4: Rate constants (ns−1) for urea dissolution in different water models at 300, 320,and 340 K. The rate constants for cubic and tablet-shaped crystals were obtained from fits tothe cube root and square root laws, respectively.Figure 3.9: The radial density profile (molecules nm−3) about the center of a cubic crystalin SPC water, measured at 0.75 ns. The dashed horizontal line in the inset indicates thesaturation density.323.4. DiscussionThe radial density profile of urea about the center of a dissolving cubic crystal is shown inFigure 3.9. We see that the urea density is nearly uniform in solution, and is much less than thesaturation density. A density gradient occurs only very near the crystal surface, suggestingthat the diffusion of urea into the bulk solution is not the rate determining step in thedissolution process. The rates we observe are qualitatively inconsistent with the observationof such a narrow diffusion layer.6We also examined the thickness of the layer that would be required to explain the observedrate. We calculated ∆N/∆t for the cubic crystal in the fixed rate law stage at ∼ 0.75 ns,and estimated dN/dt ' –130 molecules ns−1. The diffusion coefficient for urea, D ' 2.34 ×10−9 m2 s−1, and the surface area of the solid at 0.75 ns is ∼ 63.62 nm2. The saturationconcentration Cs, is 10.82 molecules nm−3. If there is a diffusion layer, the dissolution rateis given by, dN/dt = −DSCS/δ, implying that there has to be a diffusion layer of thickness∼ 12.4 nm. This is obviously too thick to be relevant in our studies, as the crystal has aradius of ∼ 2.25 nm. Thus, we can safely conclude that detachment and not solute diffusiondetermines the dissolution rate.3.4.2 Temperature dependence and activation energiesFigure 3.10: Fits to the cube root law for a cubic crystal in SPC water at 300, 320, and 340K.333.4. DiscussionWe would expect dissolution to be an activated process, and, in order to determine activationenergies, simulations for the cubic crystal were carried out at three different temperatures,300, 320, and 340 K. Rate constants (Table 3.4) were found by fitting the cube root law inthe fixed rate law regime (excluding the initial and final region of the dissolution profiles).The cube root law gives an excellent representation of the dissolution profiles, as an example,the fits for the SPC model are shown in Figure 3.10. We also used rate constants obtainedfrom the square root and the linear law to estimate the activation energy, and the energiesare listed in Table 3.5.Figure 3.11: Fit to the Arrhenius equation of the rate constants obtained for a cubic crytalin SPC water at 300, 320, and 340 K.For all water models considered, plots of of ln(k) vs. 1/T exhibit Arrhenius behavior(ln(k) = ln(A) − Ea/RT ), as illustrated in Figure 3.11 for the SPC model. The activationenergies of urea dissolution obtained for the different water models are listed in Table 3.5. Thequoted standard deviations in the activation energies were obtained by dividing each profileinto four equal parts and fitting to the rate law, such as to obtain four different estimates ofthe activation energy. We note from Table 3.5 that the activation energies associated withdissolution show no significant dependence on the water model employed, with the values forall models lying within overlapping standard deviations. For example, the activation energyobtained with the SPC/E model is nominally lower than that of SPC, yet the rate constant343.4. Discussionobtained with SPC/E at 300 K is approximately only two-thirds that of SPC. The similarityof the dissolution activation energies suggests that some other effect is influencing the pre-exponential factor in the Arrhenius equation, and giving rise to the model dependence of theactivation energies.Water Model Dissolution DiffusionCube root law Square root law linear lawTIP4P/2005 32.21 ± 3.99 – – 16.56 ± 0.30TIP3P 29.51 ± 1.91 – – 12.00 ± 0.37SPC/E 32.24 ± 3.16 – – 12.71 ± 0.43SPC 32.88 ± 3.75 32.92 ± 2.01 33.03 ± 1.05 12.15 ± 0.28Table 3.5: Activation energies (kJ mol−1) for cubic crystal dissolution and urea diffusion.Figure 3.12: Diffusion coefficients of urea (top panel) and water (bottom panel) in solution.Experimental trendlines for urea60 and water61 are included.353.4. DiscussionOne possibility we considered is that the diffusion rate of urea, or water (possibly both),is having some effect on the dissolution rate. To check this possibility, we calculated thediffusion coefficients of urea and water from mean square displacements at 300, 320, and 340K, in equilibrium solutions after the entire crystal was dissolved. The activation energies forurea diffusion estimated from Arrhenius plots (ln(D) = ln(A)−Ea/RT ) are included in Table3.5, and we note that the activation energies for urea diffusion are much smaller than thecorresponding values for dissolution. The diffusion coefficients for all four water models areplotted versus temperature in Figure 3.12, together with trendlines based on experimentaldata.60,61 We notice that for urea the TIP4P/2005 results lie closest to the experimentaltrendline, whereas for water SPC/E agrees best with experiment, with TIP4P/2005 a closesecond.From Figure 3.12, we see that for both urea and water the diffusion coefficients follow theorder, TIP3P > SPC > SPC/E > TIP4P/2005), except at 340 K where the water diffusioncoefficients for SPC/E and TIP4P/2005 are practically identical. We note that the order ofurea diffusion coefficients for the three-point water models (TIP3P > SPC > SPC/E) matchesthe variation in the dissolution rates (Figure 3.2, and Table 3.4), which suggests that differingdiffusion rates might be contributing to the different dissolution rates. However, the resultsfor TIP4P/2005 do not fit this picture. The diffusion coefficients for TIP4P/2005 are similarto those of SPC/E, but for TIP4P/2005 the dissolution rate is considerably faster at all threetemperatures. This suggests that some other, perhaps structural, effects are contributingto the different dissolution rates, but despite examining some possibilities such as radialdistribution functions, urea-water hydrogen bonding, urea-water electrostatic interactions,and urea solvation energies, we could not identify any clear correlation that would explain allof the model-dependent variation in the dissolution rates.363.5. Summary and conclusions3.5 Summary and conclusionsEmploying MD simulations we have investigated in detail the dissolution of urea nanocrystalsunder sink conditions. Both cubic and tablet-shaped crystals were considered at temperaturesof 300, 320, and 340 K. The GAFF forcefield,23 which importantly gives a good representationof crystalline urea. was used together with four common water models. The qualitativedissolution mechanism was the same for all four water models considered, but there weresurprisingly large differences in the dissolution rate.We found that all urea crystals dissolved in three stages, similar to those observed forionic NaCl nanocrystals.6 Initially, loosely bound molecules located at corners and/or edgesleave the crystal. After this process, the crystal takes on a solution annealed shape, whichpersists for most of the dissolution process during the stage we refer to as the fixed ratelaw regime. In the final stage, the dissolution profile departs from the governing rate law,crystalline structure is lost, and eventually the crystal completely vanishes into solution.During the middle stage of dissolution, the rate is well described by a classical expression,which assumes that the rate is proportional to the active surface area of the crystal. Moleculesare assumed to depart uniformly from the active surface area. After the initial stage, the activesurface area of the initially cubic crystal is essentially spherical in shape, which leads to anintegrated rate law of the cube root form. We show that the cube root law gives good fits tothe simulation results for cubic crystals, over the intermediate stage of dissolution. We alsoshow that the detachment of molecules from the crystal, and not the formation of a diffusionlayer, is the rate determining step for the dissolution process. For cubic crystals, the rateconstants obtained at different temperatures closely follow the Arrhenius equation, giving anactivation energy of ∼ 32 kJ mol−1. The activation energies obtained with all four watermolecules are very similar, indicating that the differences in rate constant come through thepre-exponential factor in the Arrhenius equation. Our investigation of this issue suggestedthat some (but definitely not all) of the model dependence observed in the rate constantsmight be connected with the model dependence of the diffusion coefficients of water and/or373.5. Summary and conclusionsurea.For the tablet-shaped nanocrystal, the overall dissolution mechanism was similar to thatof cubic crystals, but in our simulations molecules only left the tablet from the curved surface,and never from the flat faces. Thus, we identified the active surface as cylindrical in shape,which gives a square root form for the integrated rate law. The square root law gave anexcellent fit to the middle stage simulation results for the tablet-shaped crystal. However, wenote that on the time scales of our simulations, there is not much difference between cuberoot and square root rate laws, and both give good fits to the simulation profiles.Finally, we remark that the present simulations taken together with earlier work on thedissolution of NaCl crystals,6 shows that at least for some crystals dissolution might be ex-pected to follow a relatively uncomplicated three-stage mechanism. Importantly, the longmiddle stage of the dissolution process, during which the vast bulk of the crystal dissolves,is very well described by classical rate laws that depend only on the geometry of the activesurface area, as determined after the crystal has solution annealed during the initial stage.This suggests that in some cases crystal dissolution studies may not require very long simu-lations. If the mechanism we have found for urea and NaCl applies, then it would be onlynecessary to simulate the initial stage of dissolution. This would allow one to identify thegeometry of the active surface, and hence determine the rate law that governs the bulk ofthe dissolution process. More simulations of complex molecular solids are necessary in orderto assess the generality of our observations, but the fact that similar mechanisms apply tocrystals as different as NaCl and urea is an encouraging start.38Chapter 4Dissolution of Aspirin Nanocrystals4.1 OverviewThe possibility of using classical dissolution models to predict dissolution rates of molecularsolids is discussed in this chapter. Molecular dynamics simulation studies of dissolution ofaspirin nanocrystals in water are reported. All the simulations were performed in a mannerdesigned to prevent the build up of any significant concentration of aspirin in solution, andthus have a continuous dissolution of the aspirin nanocrystal. Cubic and cylindrical shapedcrystals are studied to examine the effect of shape on dissolution. The effect of temperatureis also examined. The dissolution is found to occur in four stages: an initial stage wherethe loosely bound molecules at the edges and corners go into the solution, after which therounded cubic crystal transforms to a nearly cylindrically shaped crystal in a transformationstage. The transformation stage is absent in the dissolution of cylindrical crystals. After thisstage the crystals dissolves at a nearly fixed rate, where the molecules leave the active surface(curved surface of cylindrical shape) nearly evenly, until the final dissolution stage, where thecrystal dissolves rapidly because of its increasing instability due to its reduced size. The fixedrate stage of the dissolution, which is a major fraction of the whole dissolution run, is welldescribed by classical rate equations which assume that the rate of dissolution is proportionalto the active surface area of the crystal. The detachment of the molecules from the crystal isfound to be the rate determining step for dissolution.394.2. Introduction4.2 IntroductionIn our continuing efforts to elucidate the dissolution process of molecular solids, we chooseaspirin as the model molecule for further investigations. Aspirin is an Active PharmaceuticalIngredient (API), which has been studied extensively since it’s discovery. It along with mostof the API’s are known to have low solubility in water.62 The dissolution of these solids influ-ence important processes, such as drug delivery,10 thus these are studies of high importance.However, due to poor solubility, simulation of the dissolution kinetics is very challenging.There have been few studies which have investigated the dissolution of some APIs, butsuch studies only consider distinct faces in contact with the solvent. A study of acetaminophendissolution in water by Gao and Olsen16 revealed the significance of edges and corners in theinitial stages of dissolution. So, a dissolution study of any solid would be incomplete withoutconsidering edges and corners. Elts et al.17,63 investigated the dissolution of aspirin in water,and proposed an approach that was a combination of Molecular Dynamics (MD) and kineticMonte Carlo (kMC) simulations. They used this approach to extend the time scales accessibleto simulation, and found face displacement velocities of aspirin crystals.Another recent study of the dissolution of the ionic crystal NaCl in water by Patey andLanaro6 used classical dissolution models1,2, 4, 5 to interpret dissolution rates. They show thatthese models describe the dissolution profile of the crystal to a very good extent. We took thesame approach in our study of urea in water described in Chapter 3, and found that we couldpredict the dissolution rate by considering only a fraction of the total dissolution run, byfitting the profiles to rate laws3,46,47 obtained from classical models. However before comingto a general conclusion, we test the method for a molecular solid with very low solubility inwater.In this chapter we investigate the dissolution of aspirin nanocrystals in water to furthersupport our conclusions based on the dissolution of urea nanocrystals. Aspirin was chosenbecause of it’s relatively simple structure, contrasting low solubility with respect to urea inwater, and the availability of a force field model, which replicates both solid and solution404.3. Models and methodsphase properties. The complete dissolution of aspirin nanocrystal is studied and detaileddescriptions of the various stages are presented. The effect of temperature and the initialshape of crystal on the dissolution rate is also investigated. The activation energy for thedetachment process which is the rate determining step for aspirin dissolution is estimated.The remainder of the chapter is organised as follows, models and methods used are givenin section 4.3,detailed results are described and discussed in section 4.4, and conclusions issummarized in section Models and methodsWe report various simulations of aspirin dissolution in water. In all simulations, the non-bonded interaction contains both the Lennard-Jones (LJ) and the electrostatic term as inEquation (3.1). The force-field parameters for aspirin were taken from the SwissParam64server, representing CHARMM-compatible parameters based on the Merlock Molecular ForceField (MMFF).65 These parameters produced a stable crystal lattice, and the transition ofthe solid to the solution and vice versa could be describe qualitatively. The TIP3P20 watermodel was used, as it gives the best agreement with experimental results for CHARMM modelof aspirin crystals, and has been used in previous dissolution studies of CHARMM model ofaspirin crystals.We considered various shapes (Cubic and Cylindrical) of the Aspirin crystal to investigatethe influence of shape on the dissolution rate, as used in earlier dissolution studies. The cubiccrystal was generated by repeating the Aspirin unit cell in the directions of the three axes, andthe cylindrical crystals were obtained by taking a cubic crystal and removing the moleculesoutside a circle centered about the symmetry axis of the cylindrical crystal.Cubic crystals: These contained 800 molecules of aspirin and were generated by repeat-ing the unit cell 5, 8 and 5 times along the different crystal axes. A cross section of the crystalalong the (010) crystal plane, together with other crystal planes (100 and 001), are shown inFigure 4.1. A smaller cubic crystal containing 448 molecules was also used to test our method414.3. Models and methodsfor keeping the solution undersaturated.Figure 4.1: Schmatic representation of the cross section of a cubic crystal along the (010)crstal plane. The (100) and (001) planes lie perpendicular to the plane of the paper. Theoxygen atoms are represented in red, hydrogen atoms in gray, and carbon atoms in black,respectively.Cylindrical crystals: We investigated three different cylindrical crystals with differentflat surfaces (the (100), (010) and (001) faces of the aspirin crystal). These contained 550,464 and 490 molecules, respectively. Cylindrical crystals were generated by taking the cubiccrystal and considering the three crystal axes as the axis for the three different cylindricalcrystals, and then removing the molecules outside a circle of radius 2.4 nm centered about thesymmetry axis of the cylindrical crystal. Top views (flat surface) and the side views (curvedsurface) of the three crystals are shown in Figure 4.2.424.3. Models and methodsFigure 4.2: a) Side views (curved surface), and b) Top views (flat surface) of the cylindricalcrystals. The crystals with (001), (010) and (100) faces as the flat surface are shown in i), ii)and iii), respectively. The oxygen atoms are represented in red, hydrogen atoms in gray, andcarbon atoms in black, respectively.All the simulations were carried out under NPT conditions with the pressure fixed at 1bar. The details for the different initial configuration (temperature of the system and thenumber of aspirin and water molecules) in the various simulations are summarized in Table4.1.434.3. Models and methodsUrea Water Temp(K)Cube 800 39688 300800 39688 320800 39688 340Cylinder 464 37970 300490 37445 300550 37442 300Table 4.1: A summary of the number of aspirin and water molecules used, shapes of the initialcrystals, and the temperatures at which simulations were performedAll simulations were performed using the GROMACS molecular dynamics package,59 ver-sion package 4.5.5. The temperature was controlled using the Nose´-Hoover thermostat35,36with a time constant of 2.0 ps. The pressure was maintained by using the Parrinello-Rahmanbaraostat37 with a compressibility of 4.46 × 10−5 bar−1 and a time constant of 2.0 ps. Thetime step for all the simulations was 1 fs. Periodic boundary conditions were used and theshort range interactions were truncated at 1 nm. The electrostatic interactions were takeninto account using the fourth order Particle-Mesh Ewald method33 with a Fourier spacing of0.12nm. All the bonds were constrained using the LINCS algorithm.38The initial configurations for the simulations were generated by placing the crystal in thecenter of a cubic box of length 11 nm, and the box was filled with water molecules. The watermolecules within the crystal lattice were then removed, and the system was relaxed for 0.2 nskeeping the aspirin molecules fixed, in order to equilibriate the solvent molecules around thecrystal.The experimental saturation concentration66 of aspirin in water at 300-340 K correspondto aspirin mole fractions in the range 0.00112 - 0.0042. The extremely low solubility of as-pirin results in the solution becoming oversaturated after just a small fraction of the crystaldissolved, and a new procedure was required to keep the solution undersatured. Earlier stud-ies17,63 used a slab of sticky dummy atoms to irreversibly trap molecules in solution (within acutoff radius of the dummy atoms), by using a potential well. The trapped molecules were re-moved permanently from the system and constant undersaturation was observed throughoutthe dissolution. We devised another simple method to ensure undersaturation at all instances444.3. Models and methodsduring the simulation, resulting in continuous dissolution of the aspirin crystal.Aspirin molecules which have passed into solution from the crystal and are outside asphere (centered at the center of mass of the crystal) of fixed radius are removed at regularintervals. The box was then translated to make the center of the box coincide with the centerof the aspirin crystal. This was done to ensure that the sphere always remained inside thebox to prevent removal of any misrepresented aspirin molecules due to the periodic boundaryconditions. The system was then allowed to evolve for another time interval, and the wholeprocess for removal of aspirin molecules was repeated until complete dissolution of the crystalwas achieved. The number of aspirin molecules leaving the crystal during different timeintervals was monitored, and a time interval which ensured undersaturation at all instancesduring the simulation was chosen. The time interval used for removal of aspirin molecules was2 ns for simulations at 300 K, and 1 ns for simulations at 320 and 340 K. The actual valueof the radius did not greatly affect our results, as the dissolution profiles were qualitativelysimilar for values of 3.5 and 4.5 nm, as shown in Figure 4.3. The minor deviation at the end ofthe dissolution profile is mostly stochastic in nature, since very similar profiles were obtainedwhen the radii were reversed at 240 ns (see inset in Figure 4.3). In our simulations we chosea value of 4.5 nm for the radius of the sphere beyond which aspirin molecules in solution wereremoved. The simulations required from 100 to 620 ns for complete dissolution, dependingon the shape of the crystal and the temperature of the system.The molecules were identified as a part of the crystal or the solution by means of a simpleorder parameter based on the number of neighbors. The number of neighbors was identifiedby counting the number of aspirin molecules within a sphere of radius 1.10 nm centered at thecenter of mass of the molecule. Molecules that were in the crystal had neighbors in the range15-22, and the molecules in the solution had 1-3 neighbors. So, we chose the intermediatenumber 8 to be our cutoff value. A molecule with 7 or fewer neighbors was considered to bea part of the solution, and any molecule with 8 or more neighbors was identified as part ofthe crystal. The radius of the sphere for counting the neighbors and the cutoff value did notaffect our results significantly, as the profiles were qualitatively similar for values of 0.70 and454.4. Results and discussion1.40 nm.Figure 4.3: Dissolution profiles for the cubic crystal with different radii of the sphere used forremoval of molecules in the solution. The profiles shown in the inset are for the same crystalswith the sphere radii reversed at 240 ns in the dissolution runs.4.4 Results and discussion4.4.1 Stages of dissolutionThe dissolution is found to occur in four different stages. The initial stage can be dividedinto two sub-stages, a very rapid dissolution of the loosely bound molecules, and a slowerdissolution during which the crystal take a solution annealed shape. This is followed by anearly steady dissolution rate for most of the dissolution run, after which the crystal dissolvesat a very rapid rate as it becomes increasingly unstable due to its decreasing size.A. Initial stageDuring this phase of the dissolution, the loosely bound molecules at the corners and edgesof the crystal get detached and passes into solution. This process depends on the shape of464.4. Results and discussionthe crystal, as the shape determines the number of molecules on the edges and corners of thecrystal, and leads to exposure of more crystal molecules to the solution, thus enabling thesolution to anneal the crystal to a shape that gives an active surface from which the moleculesleave at a nearly even rate.Figure 4.4: The dissolution profile for a cubic crystal displayed as the number of molecules vstime. The profiles for a cubic crystal and a cylindrical crystal with the (100) face as its flatsurface are shown in the inset.A dissolution profile for an initially cubic crystal is plotted in Figure 4.4, and snapshotsof the crystal at various points in the dissolution run are shown in Figure 4.5. For the cubiccrystal, as can be seen from the dissolution profile in Figure 4.4, the crystal loses nearly 100molecules in the first 20 ns of the dissolution run. It is also evident from the snapshots inFigure 4.5 that during the first 20 ns the crystal has lost nearly all of the molecules at theedges and corners. After these initial detachments, the crystal still has a basically cubic shapewith six faces and rounded edges and corners. It takes the cubic crystal another 100 ns totake a cylindrical shape, which is the next phase of the dissolution, and is detailed in the nextsection.474.4. Results and discussionFigure 4.5: Snapshots of the (100) face corresponding to different points in the dissolutionprofile of a cubic crystal. Each red sphere indicates an aspirin molecule.Figure 4.6: Dissolution profiles for three cylindrical crystals with different flat surfaces.484.4. Results and discussionFor the three different cylindrical crystals, there is a similar trend in the dissolution profile,where they lose the molecules on the edges of the crystal very rapidly in the first part of thedissolution run. This is evident from the dissolution profiles in Figure 4.6. Also it can be seenfrom the snapshots of the cylindrical crystal in Figure 4.7, that the crystal loses moleculesfrom the curved surface and not the flat face. This stage is followed by a nearly fixed ratedissolution, where molecules come off the crystal evenly from the active surface area.Figure 4.7: Snapshots corresponding to different points in the dissolution profile of a cylin-drical crystal with the (100) face as its flat surface. Each red sphere indicates an aspirinmolecule.B. Transformation stageThe rounded cubic crystal after the initial rapid dissolution stage still has the three crystalfaces (100, 001, and 010) that were present in the initial crystal. These surfaces have differentdisplacement velocities,17 with the (100) face having the smallest, and the other faces (010)and (001) having similar velocities. This makes the (100) face the most stable surface of allthe surfaces exposed to solution. The water gradually consumes molecules from corners of theother surfaces at a very rapid rate, in comparison with molecules from the (100) face. Thisresults in the formation of a nearly cylindrical shaped crystal with the (100) face as the flatsurface, and the other faces being less stable, are attacked to form a nearly uniform curvedsurface. This transformation of the rounded cubic crystal to a nearly cylindrical crystal takesaround 100 ns, and can be seen from the snapshots of the dissolution run shown in Figure4.5. Once the cylindrical crystal is formed in solution, molecules start to come off the curvedsurface evenly resulting in a nearly fixed rate of dissolution, and this cylindrical shape of the494.4. Results and discussioncrystal is retained for most of the dissolution run, as evident from Figure 4.5.The active surface for cylindrical crystals are the curved surfaces from the beginning ofthe dissolution run, which is the solution annealed stable shape for these crystals. This isthe reason we do not see a transformation stage in the cylindrical crystals, and they directlytransition to the fixed rate law stage after the initial stage.C. Fixed rate law stageAfter the initial stages of dissolution, all the crystals follow a similar trend. The dissolutionprofile remains almost linear during this stage, which is a major part of the whole dissolutionrun. The crystal retains the stable shape (cylindrical), attained after the initial stages of thedissolution run, throughout this stage.During this phase of dissolution, molecules leave the crystal evenly from the active surfacearea formed after the initial stages. The curved surface of the cylindrical crystals is iden-tified as the active surface in these crystal, as no molecules leave the crystal from the flatsurfaces. This essentially uniform detachment of molecules from the active surface results ina dissolution process which follows a fixed rate law, until the final rapid dissolution stage.Next we evaluate the effect of crystal shapes on the dissolution rate and also compare thedissolution profile in the fixed rate law stage with different classical models mentioned above.D. Influence of crystal shapeCubic and cylindrically shaped crystals were used to investigate the effects of shape. It isevident from the dissolution profiles (Figures 4.4 and 4.6) that all the crystals have similarinitial and final stages, where there is a rapid dissolution of the molecules from the crystal(details in sections 4.4.1.A and 4.4.1.F, above). Cylindrical crystal do not undergo the trans-formation stage, where cubic crystals transform from a rounded cubic structure to a nearlycylindrical structure (see section 4.4.1.B for details).Once this transformation is complete, all crystals have a nearly fixed rate of dissolution.It can be seen from the inset in Figure 4.4 that the dissolution profile for the fixed rate stage504.4. Results and discussionof a cylindrical crystal with the (100) face as the flat surface is qualitatively similar to that ofa cubic crystal. An interesting observation from these profiles is that the rate of dissolutionof a crystal with more than one flat surface (having different displacement velocities) is verysimilar to a cylindrical crystal with the most stable face of the crystal as the flat surfaces.E. Comparison with dissolution modelsAs discussed before, there have been a number of dissolution studies over the last century,and several models in the form of differential and integral rate laws have been proposed. Weemploy a very simple rate law to describe the fixed rate stage of these dissolution profiles.Since all of our simulations have been carried out under sink conditions, the detachmentprocess is expected to be the rate determining step, and the rate can be described asdN(t)dt= −kSactive(t), (4.1)where N(t) is the number of molecules in the crystal at time t, Sactive(t) is the active surfacearea where the detachment process takes place, and k is a constant. This law, which takesinto account the changing surface area, was first suggested by Brunner and Tolloczko.2The differential form is used to get the integral rate law using different forms of theactive surface area, Sactive(t). The active surface area in our case is the curved surface ofthe cylindrical structure. Thus, using area to volume relationships gives us the relationSactive(t) ∝ N1/2(t), which leads to the integrated rate law√N(t) =√N0 − kt, (4.2)where k is the rate constant. The fixed rate regimes for the crystals were then compared withthe law, by fitting the integrated rate law equation to the dissolution profiles. The goodnessof fit parameters, R2, for each crystal are listed in Table 4.2. Results for the cube root andthe linear rate laws are also included in the table. The best fits clearly corresponds to the ratelaw derived above for the cylindrical crystals, while in case of the cubic crystal the differences514.4. Results and discussionin the goodness of fit parameters are too small for us to draw any firm conclusion.Cubic Cylindrical(100) face (010) face (001) facelinear 0.9892 0.9922 0.9940 0.9937square root 0.9852 0.9950 0.9941 0.9976cube root 0.9829 0.9947 0.9930 0.9976Table 4.2: Goodness of fit parameters R2 obtained for different rate lawsWe also argue that the diffusion of the molecules in the solution is not the rate determiningstep, as the simulation was done in a way to prevent build up of any diffusion layer aroundthe crystal. We have calculated the diffusion coefficients of aspirin in the solution at 300, 320,and 340 K, and the values are listed in Table 4.3. The activation energy for aspirin diffusionwas estimated from the Arrhenius equation ( ln(D) = ln(A)− Ea/RT ), and is 27.82± 4.24kJ mol−1. It can be noted that the value is significantly smaller than the activation energyfor dissolution (48.31 ± 6.63 kJ mol−1, see section 4.4.2, below). Thus, the diffusion layermodel does not apply to the dissolution of the Aspirin.Diffusion Coeffecients (cm2 s−1× 10 −5)TIP3P TIP4P/2005300 K 0.88 ± 0.53 0.78 ± 0.25320 K 1.90 ± 0.50 1.43 ± 0.43340 K 3.24 ± 0.67 1.97 ± 0.85Table 4.3: Diffusion coefficients of aspirin in solution.Figure 4.8: Snapshots of the last stages in the dissolution profile of a cubic crystal. Each redsphere indicates an aspirin molecule.524.4. Results and discussionF. Final stageDissolution profiles in Figure 4.4, suggest that there is a sharp increase in the dissolutionrate during the final stages of the dissolution run. The crystal becomes very unstable, anddisappears rapidly once it has been reduced to a very small size. During this phase the crystalstarts to loose its crystalline structure and then disintegrates completely. The is evident inthe snapshots shown in Figure Temperature dependence and activation energyThe detachment of the molecules from the crystal is expected to be an activated process,and to estimate the activation energy we carried out simulations of the cubic crystal at threedifferent temperatures, 300, 320, and 340 K. The rate constants were calculated by fitting thefixed rate stage of the dissolution profiles to the square root law. The profiles fit the rate lawexcellently, as shown in Figure 4.9. The rates follow the Arrhenius equation (Figure 4.10)ln(k) = ln(A)− Ea/RT , (4.3)giving an activation enrgy of 48.31 ± 6.63 kJ mol−1.Figure 4.9: Fits to the square root law for a cubic crystal at temperatures of 300, 320, and340 K.534.4. Results and discussionFigure 4.10: Fits to the Arrhenius equation of the rate constants obtained for a cubic crystalat temperatures of 300, 320, and 340 K.The standard deviations shown as error bars in Figure 4.10 were calculated by dividingthe fixed rate stages of the dissolution profiles into five equal parts, and rate constants werecalculated by fitting each part to the square root law. These five different rate constant werethen used to calculate the standard deviation, shown as error bars in Figure 4.10. Similarly,the five different sets of rate constants were then used to calculate the standard deviation inactivation energy.4.4.3 Model dependenceThe dissolution of aspirin nanocrystals in water was also studied using the TIP4P/200521water model. The qualitative dissolution profile for the cubic crystal does not depend on thewater model used, and we see the four stage dissolution in both cases. However, the dissolutionrate is found to depend on the water model used, with the dissolution in the TIP4P/2005water model occuring considerably slower than in TIP3P water. The activation energy fordissolution in TIP4P/2005 water model was estimated using the Arrhenius behavior of therate constants at 320 and 340 K, and was found to be 52.72 kJ mol−1. The plot of thedissolution profile in TIP4P/2005 water model is shown in Figure 4.11, and it can be noted544.5. Summary and conclusionsthat the profiles are well fit by the square root law.Figure 4.11: Fits to the square root law for a cubic crystal in TIP4P/2005 water at temper-atures of 300, 320, and 340 K.The dependence of the dissolution rate on the water model used can be explained by theslight difference in activation energy for the water models, despite the high standard deviationin the estimation of the activation energy. We assume that the activation energy is dependenton both the crystal and crystal-water interface, though a very weak dependence on the latter.This is justified as we found that the aspirin molecules interact with TIP4P/2005 water modelmore strongly, thus they may form a more stable crystal-water interface in TIP4P/2005 watermodel, leading to a slightly bigger activation barrier, and a slower dissolution rate. Otherfactors such as solvation time, aspirin-water hydrogen bonding, and diffusion coefficients (seeTable 4.3), are very similar for both the cases.4.5 Summary and conclusionsMolecular dynamics simulations were carried out to study the dissolution of a molecular modelfor aspirin nanocrystals of different shapes and size in detail. Cubic and cylindrically shapedcrystals were used in simulations, where all dissolution was performed under sink conditions.In our studies, the dissolution was found to occur in four stages for cubic crystals and three554.5. Summary and conclusionsstages for cylindrical crystals. The initial stage, where loosely bound molecules at the cornersand edges of the crystal pass into the solution, was common in all the simulations. Thisstage depends on the initial structure of the crystal as the structure influences the number ofmolecules on the edges/corners of the crystal. After the initial stage, both crystals retainedtheir initial shape with the edges being smoothed by the solution. A transformation stage,which occurred only in the dissolution of cubic crystals, transformed the well rounded cubiccrystal into a nearly cylindrical shape. Once both crystals had acquired a nearly cylindricalshape, molecules left the crystal from the curved surface of the cylindrical crystal essentiallyevenly, giving a nearly fixed dissolution rate until the final stage of dissolution, where thecrystal dissolves at a very rapid rate because of its increasing instability due to the reductionin size. The stages appears to be similar to that of the ionic NaCl nanocrystals and molecularurea nanocrystals.The fixed rate stage of the dissolution was well described by considering the rate to beproportional to the active surface area of the crystal. For both the crystals, the active surfacearea was the curved surface of the cylindrical shape they attain after the initial stage. Thuswe obtain the square root law, which is appropriate for crystals having cylindrical shape.We show that the square root law describes the fixed rate stage very well, by examining thegoodness of the fit parameter for each simulation. We also show that the detachment ofmolecules from the crystal is the rate determining step, and not the formation of a diffusionlayer around the crystal. We determined the activation energy for the detachment process byfitting the rates of dissolution of the cubic crystal at three different temperatures 300, 320 and340 k to the Arrhenius equation. The activation energy was approximately 48.31 kJ mol-1,which is significantly bigger than the value obtained for urea dissolution (∼ 32 kJ mol−1).We conclude from our simulations of urea and aspirin nanocrystals, that the classicalmodels give an accurate description of the fixed rate stage of dissolution of these nanocrystal.Thus we can use these models to study the dissolution of different molecular nanocrystals,and predict the rates by running only a fraction of the simulation, therefore reducing theproblem of the long time scales required for a dissolution run.56Chapter 5Summary and Conclusions5.1 SummaryDissolution research has been developing for over a century as a field in physical chemistry, anddifferent models have been put forward as a description of the dissolution process. However,in the last few decades dissolution of microscopic molecular crystals has gained considerableinterest, because of the importance of dissolution properties in drug bioavailability. How-ever, most of the active pharmaceutical ingredients exhibit poor solubility or permeability,and thus simulations of the dissolution process is very challenging, even after the immenseincrease in computational power. In this thesis, we investigated the dissolution process ofmolecular solids, using urea and aspirin nanocrystals as models, and propose a method topredict dissolution rates by analyzing the results of the simulations.In Chapter 2, the models used to represent urea, aspirin and water molecules in thesimulation are presented, as well as the different algorithms and methods of the moleculardynamics simulations.Chapter 3 describes an investigation of the dissolution process of urea nanocrystals inwater. We found that the dissolution of urea can be described as a three step process, wherethe crystal first looses the edges and corners very rapidly, followed by a steady dissolutionrate until it reaches a certain size (. 200 molecules), and then losses its crystalline structureand dissolves completely into solution. It was also found that the rate laws obtained fromthe classical models described the fixed rate regime of the dissolution profile to a very goodextent. Additionally, dissolution was found to be an activated process, with detachment ofmolecules from the crystal being the rate determining step. The activation energy for urea575.2. Future directionsdissolution was found to be ∼ 32 kJ mol−1 and was similar for different water models used.However, the dissolution rate varied across the water models, implying some dependence onthe pre-exponential factor in the Arrhenius equation.In Chapter 4 we further investigated the dissolution process by using aspirin nanocrystalsas a model. The dissolution was found to occur in four steps, with an additional transformationstage where the crystal, with molecules at the edges and corners gone into the solution,dissolves slowly to take a solution annealed shape before the steady dissolution regime. Therate law obtained from the classical model was found to describe the steady part of thedissolution profile again to a very good approximation. The activation energy was found tobe ∼ 49 kJ mol−1, which is perhaps expected given the low solubility of aspirin as comparedwith urea in water.Based on our simulation results, it can be concluded that classical rate laws can be used topredict dissolution rates for some molecular solids, by fitting an early part of the dissolutionprofile to simple rate law expressions.5.2 Future directionsOur results show a surprisingly strong dependence of the dissolution rate on the water modelused for simulation. This implies that small differences in the electrostatic charges on the watermodel influences the dissolution rate, and since the majority of the dissolution processes ofinterest involve complex mixtures of different solvents (water, acids, and bases), it would beof interest to examine if the proposed dissolution mechanism is valid in these more complexsituations.A partial explanation for the difference in dissolution rates came from the difference indiffusion coefficients. We found that the diffusion coefficients of urea, and water models aredirectly proportional to urea dissolution rates in three-point water models. However, this didnot explain the dissolution rates in the four-point water model. Another explanation was thedifference in activation energies, which could lead to different rates, but the large standard585.2. Future directionsdeviations in the activation energies makes it difficult to evaluate this possible explanation.Thus, more studies are definitely required to study this difference in dissolution rate acrossmodels.59Bibliography[1] A. A. Noyes and W. R. Whitney, J. Am. Chem. Soc. 19, 930 (1897).[2] L. Bruner and S. Tolloczko, Z. Phys. Chem. 283 (1900).[3] A. Hixson and J. Crowell, Ind. Eng. Chem. 23, 1160 (1931).[4] E. Brunner, Z. phys. Chem. 43, 56 (1904).[5] W. Nernst, Z. phys. Chem. 47, 52 (1904).[6] G. Lanaro and G. N. Patey, J. Phys. Chem. B 119, 4275 (2015).[7] M. Wilderman, Z. Phys. Chem. 66, 445 (1909).[8] A. B. Zdanovskii, Zhur. Fiz. Khim. 20, 869(1946).[9] P. V. Danckwerts, Ind. Eng. Chem. 43, 1460 (1951).[10] A. Dokoumetzidis, and P. Macheras, Int. J. Pharma. 321, 1 (2006).[11] S. S. Jambhekar and P. J. Breen, Drug Discovery Today 18, 1173 (2013).[12] J. B. Dressman, G. L. Amidon, C. Reppas and V. P. Shah, Pharm. Res. 15, 11 (1998).[13] O. Anand, X. Y. Lawrence, D. P. Conner and B. M. Davit, AAPS J. 13, 328 (2011).[14] J. Siepmann and F. Siepmann, Int. J. Pharma. 453, 12 (2013).[15] S. Piana and J. D.Gale, J. Am. Chem. Soc. 127, 1975 (2005).[16] Y. Gao and K. W. Olsen, Mol. Pharmaceutics 10, 905 (2013).60Bibliography[17] E. Elts, M. Greiner and H. Briesen, Crys. Growth Des. 7, 4154 (2016).[18] H. J. Berendsen, J. P. Postma, W. F. van Gunsteren and J. Hermans, Intermolecularforces, Springer:Dordrecht, Netherlands, 331 (1981).[19] H. J. Berendsen, J. R. Grigera and T. P. Straatsma, J. Phys. Chem. 91, 6269 (1987).[20] W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey and M. L. Klein, J.Chem. Phys. 79, 926 (1983).[21] J. L. Abascal and C. Vega, J. Chem. Phys. 123, 234505 (2005).[22] P. Xiu, Z. Yang, B. Zhou, P. Das, H. Fang and R. Zhou, J. Phys. Chem. B 115, 2988(2011).[23] G. A. O¨zpınar, W. Peukert and T. Clark, J. Mol. Model. 16, 1427 (2010).[24] G. A. O¨zpınar, F. R. Beierlein, W. Peukert, D. Zahn and T. Clark, J. Mol. Model. 18,3455 (2012).[25] P. Vaughan and J. Donohue, Acta. Crystallogr. 5, 530 (1952).[26] J. E. Worsham, H. A. Levy and S. E. Peterson, Acta. Crystallogr. 10, 319 (1957).[27] S. Swaminathan and B. M. Craven, Acta. Crystallogr. Sect. B 40 , 300 (1984).[28] C. C. Wilson, New J. Chem. 26, 1733 (2002).[29] M. Greiner, E. Elts, J. Schneider, K. Reuter and H. Briesen, J. Cryst. Growth 405,122 (2014).[30] R. W. Hockney, S. P. Goel and J. Eastwood, J. Comp. Phys. 14, 148 (1974).[31] M. P. Allen and D. J. Tildesley, Computer simulation of liquids, Oxford UniversityPress:New York, 1989.[32] P. P. Ewald, Ann. Phys. 369, 253 (1921).61Bibliography[33] T. Darden, D. York and L. Pedersen, J. Chem. Phys. 98, 10089 (1993).[34] J. W. Cooley and J. W. Tukey, Math. Comp. 19, 297 (1965).[35] S. Nose´, J. Chem. Phys. 81, 511 (1984).[36] W. G. Hoover, Phys. Rev. A 31, 1695 (1985).[37] M. Parrinello and A. Rahman, J. Appl. Phys. 52, 7182 (1981).[38] B. Hess, H. Bekker, H. J. Berendsen and J. G. Fraaije, J. Comput. Chem. 18, 1463(1997).[39] P. H. Stahl and C. G. Wermuth, Handbook of Pharmaceutical Salts: Properties, Selec-tion and Use, Verlag Helvetica Chimica Acta, Zu¨rich, 2002.[40] J. D. Wuest, Nat. Chem. 4, 74 (2012).[41] H. Crawford and J. McIntosh, Arch. Int. Med. 36, 530 (1925).[42] J. H. Meessen and H. Petersen, Ullmann’s Encyclopedia of Industrial Chemistry. WileyOnline Library, (1996).[43] M. Dunky, Int. J. Adhes. Adhes. 18, 95 (1998).[44] W. Raab, Cosmet. Toiletries 105, 97 (1990).[45] P. Costa and J. M. S. Lobo, Eur. J. Pharm. Sci. 13, 123 (2001).[46] P. Niebergall, G. Milosovich and J. Goyan, J. Pharm. Sci. 52, 236 (1963).[47] W. Higuchi and E. Hiestand, J. Pharm. Sci. 52, 67 (1963).[48] F.-M. Lee and L. E. Lahti, J. Chem. Eng. Data 17, 304 (1972).[49] R. A. Kuharski and P. J. Rossky, J. Am. Chem. Soc. 106, 5786 (1984).[50] M. Salvalaglio, T. Vetter, F. Giberti, M. Mazzotti and M. Parrinello, J. Am. Chem.Soc. 134, 17221 (2012).62Bibliography[51] J. F. Brandts and L. Hunt, J. Am. Chem. Soc. 89, 4826 (1967).[52] D. B. Watlaufer, S. K. Malik, L. Stoller and R. L. Coffin, J. Am. Chem. Soc. 86, 508(1964).[53] M J. Schick, J. Phys. Chem. 68, 3585 (1964).[54] E. M. Duffy, D. L. Severance and W. L. Jorgensen, Isr. J. Chem. 33, 323 (1993).[55] L. J. Smith,H. J. Berendsen and W. F. van Gunsteren, J. Phys. Chem. B 108, 1065(2004).[56] W. L. Jorgensen and J. D. Madura, Mol. Phys. 56, 1381 (1985).[57] C. Vega and J. L. Abascal, Phys. Chem. Chem. Phys. 13, 19663 (2011).[58] C. M. Riley and T. W. Rosanske, Development and Validation of Analytical Methods.;Elsevier Science: New York, (1996).[59] S. Pronk, S. Pa´ll, R. Schulz, P. Larsson, P. Bjelkmar, R. Apostolov, M. R. Shirts,J. C. Smith, P. M. Kasson and D. van der Spoel, Bioinformatics 29, 845 (2013).[60] A. Sadeghi, D. Kissel and M. Cabrera, Soil Sci. Soc. Am. J. 52, 46 (1988).[61] M. Holz, S. R. Heil and A. Sacco, Phys. Chem. Chem. Phys. 2, 4740 (2000).[62] N. Variankaval, A. S. Cote and M. F. Doherty, AIChE J. 54, 1682 (2008).[63] M. Greiner, E. Elts and H. Briesen, Mol. Pharm. 11, 3009 (2014).[64] V. Zoete, M. A. Cuendet, A. Grosdidier, and O. Michielin, J. Comput. Chem. 32,2359(2011).[65] T. A. Halgren, J. Comput. Chem. 17, 490(1996).[66] A. Apelblat and E. Manzurola, J. Chem. Thermodyn. 31, 85(1999).63


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