A Computational Route to Hansen Solubility Parameters of Nitrated

A COMPUTATIONAL ROUTE TO HANSEN SOLUBILITY PARAMETERS OF NITRATED SPECIES THROUGH A QUANTITATIVE STRUCTURE ACTIVITY RELATIONSHIP (QSAR) 1 CHAPTER I INTRODUCTION For many years scientists have used the “like-dissolves-like” criteria to select solvents used in their laboratories. This method is often based on trial and error and a subjective evaluation of a solvent’s polarity. This method has had varying degrees of success. One of two outcomes is usually observed. If the solvent selected provided too much solubility, then the solute becomes trapped in the solvent. This is a problem if the solute is to be extracted from the solvent later. If the solvent selected does not provide enough solubility, then the mixture is left in a heterogeneous state. This requires the addition of more solvent to achieve complete solubility. Either problem is solved with a waste of materials, particularly the solvent. In recent years scientists have endeavored to remove the guesswork from solvent selection. They have turned to model-based approaches which are capable of quantifying the solubility or solvency of a solute or solvent, and the components of their models are termed solubility parameters. History of Solubility Parameters Solubility parameters have been around for a many years and are of use to a wide variety of industries. There are many models for solubility parameters, but the way in which all solubility parameters are used is this: materials with sufficiently similar solubility parameters will interact strongly with each other. Two liquids with similar solubility parameters will be miscible. A solvent and a solute with similar solubility parameters will likely result in a mixture that is homogeneous. A liquid on a smooth, solid surface with similar solubility parameters will spread evenly and with minimal contact angle. In this way solubility parameters can provide a systematic way of 2 selecting a single solvent or a blend of solvents for the dissolution of a given solute. In the paint and polymer industry this is applicable to solvent and pigment/binder formulations. Solubility parameters have also been useful to the cleaning industry where the cleaning solution is targeted to the particular soil and surface. The term solubility parameter was first coined by Hildebrand and Scott (1,2). Hildebrand’s parameter is defined as the square root of the cohesion energy per unit of volume: ⁄ (1) Where Ec is the energy of vaporization and Vm is the molar volume. The energy of vaporization is proportional to the amount of cohesive energy of the material in the following way. The cohesive energy of a system is the sum of all intermolecular forces together. The energy required for any given molecule to escape the condensed phase and become a vapor would necessitate the breaking of all intermolecular associations. Thermodynamics require that a spontaneous process will be characterized by a negative free-energy change. As such, the thermodynamics of mixing are governed by the free-energy equation of thermodynamics: ∆ ∆ ∆ (2) Where GM is the free-energy of mixing, HM is the enthalpy of mixing, T is the absolute temperature of the system, and SM is the entropy of mixing. Hildebrand and Scott proposed an equation for the determination of the enthalpy of mixing: ∆ (3) Where 1 and 2 are the volume fractions of mixture components 1 and 2, 1 and 2 are the Hildebrand solubility parameters of components 1 and 2, and VT is the total 3 volume of the mixture. Equations 2 and 3 illustrate the fundamental usage of solubility parameters. As the Hildebrand parameters of components 1 and 2 approach the same value the enthalpy of mixing will approach zero independent of their respective volume fractions. When the solubility parameters of components 1 and 2 are sufficiently similar the enthalpy of mixing will reach a value where it can no longer dominate the entropy term of equation 2. This results in a negative free-energy of mixing and spontaneous mixing will occur. One of the fundamental weaknesses of Hildebrand’s work is that it is only applicable to “regular solutions.” Regular solutions have no “orienting or chemical effects” and all molecules are oriented and distributed randomly throughout (3). Unfortunately most chemical mixtures have dipole or hydrogen-bond orientations. Multi-component solubility parameter systems offer some measure of improvement by accounting for arrangement and orientation effects. Hansen’s widely popular three-component model of solubility parameters treats intermolecular attraction forces independently by giving each attractive force its own solubility parameter. Hansen chose three parameters for his model and gave them the symbols D, P, and H to account for dispersion (induced dipole-induced dipole) forces, polar (dipole-dipole) forces, and hydrogen-bonding forces respectively. They may be related to Hildebrand’s “total solubility parameter” by the square root of the sum of the squares of the three Hansen parameters (4). (4) It now may be easy to see in Table 1 how Hildebrand’s parameter predicts two substances such as ethylene carbonate and methanol (T = 29.6) to be miscible when they 4 are not. A glance at the Hansen solubility parameters (HSPs) of each species provides an answer. Although they possess the same total quantity of cohesive energy per unit volume, each species has that energy partitioned into different attractive forces. The fundamental usage of solubility parameters still applies, however, for two species to interact strongly they must now have three sufficiently similar parameters instead of one. Table 1. A comparison of two immiscible liquids predicted to be miscible by Hildebrand's equation. Species ethylene carbonate methanol D / MPa1/2 P / MPa1/2 H / MPa1/2 T / MPa1/2 19.4 15.1 21.7 12.3 5.1 22.3 29.6 29.6 One may think of this another way. If a species’ Hansen solubility parameters were plotted as a small sphere in a Cartesian coordinate space with its center defined as a point located by the three parameters, those species which share a significant volume within that coordinate space should be mutually soluble. Two species may reside at the same distance from the origin (represented by T) yet occupy two separate volumes of space and therefore would not be soluble. There are other multi-component solubility parameter models with varying strengths and weaknesses. All have achieved some measure of success in academia or the private sector. For example, Karger, Snyder, and Eon have developed a solubility parameter system for chromatographic solvent selection (5). Another system developed by Beerbower, Martin, and Wu provides for a more sophisticated treatment of Hansen’s hydrogen-bonding parameter (6,7). Because of this increased ability to handle the complex nature of hydrogen-bonding interactions Beerbower, Martin, and Wu’s four- 5 component system is rapidly gaining acceptance in many fields of study. However, with four components, one loses the ability to visualize the interaction in a three-dimensional Cartesian space. Previous Methods of Calculating Solubility Parameters One weakness of this model of solubility is the confusing and tedious methods used to calculate the parameters. These methods attempt to be universal but still do not apply to many cases. This is dependant upon how much information is known about the species of interest. For some species the methods yield wholly inaccurate results. It is this weakness that the work presented in this document aims to resolve. Perhaps the most straightforward way in which to determine solubility parameters is to use solvents with well determined solubility parameters to experimentally determine the parameters of the unknown. This is akin to using radio tower triangulation to locate a signal between them. Picture the Cartesian coordinate system mentioned earlier. An experiment is performed where ten solvents with well known solubility parameters are used to solvate a species of unknown parameters. The solubility of the unknown species can then be ranked by percent mass composition. The user then determines an operational definition of “good” solubility. The solute is represented as a point in space with a radius of interaction which encompasses the “good” solvents and excludes the “bad” solvents (8). This produces the most reliable parameters as they are dependant upon empirical evidence alone. The obvious drawback to this method is the impracticality and cost of applying it to a large number of species. Solute dimerization can give unusual results on rare occasions. The other obvious drawback is the 6 fundamental criteria for the experiment: the unknown species must be soluble in something. A cost-effective estimation of solubility parameters can be achieved via group contribution methods. This method operates on the principle that each functional group or part of a molecule accounts for some piece of the solubility parameters of the whole molecule. The molecule can then be built in pieces and the contributions of the individual pieces can then be summed to give the solubility parameters of the whole. This method was pioneered by van Krevelen and Hoftyzer(9) to build up the solubility parameters of polymers. It was later refined by Stefanis and Panayioutou(10). The universality of the group contribution method is often a problem. Most functional groups are not isolated from others and as such do not exist in an independent fashion. On the contrary, groups are absolutely affected by their nearest neighbors either electrostatically or through the conjugative effects of molecular orbitals. This has an effect on the overall solubility parameters of the molecule. The other phenomenon which this method of calculation does not address is the orientation of groups. This is important to the polar parameter. If two moieties have separate dipole moments but are oriented in opposite directions then they will destructively interfere. Conversely, if they are oriented in the same direction they will constructively interfere. If the two moieties provide the only significant source of dipole-dipole interaction then this method will likely produce inaccurate results for the polar solubility parameter. The use of planes of symmetry to correct for these effects introduces more complexity for very little benefit. Of the two problems, electronic and orientation, Stefanis and Panayioutou have successfully addressed the former. They have used conjugation theory to provide a more 7 sophisticated group contribution table with first order and second order groups. Unfortunately what results is a longer and more confusing process of group contribution, albeit more successful than that of van Krevelen and Hoftyzer for many molecules. The fundamental problem of group contribution methods is that by their very nature they cannot be entirely inclusive. There are an infinite number of ways in which functional groups can be assembled and this results in an infinite number of group contributions. There are ways of calculating individual solubility parameters. When possible, the dispersion parameter is calculated from the critical temperature and molar volume of an alkane homomorph in a procedure devised by Blanks and Prausnitz (11). The homomorph structure is defined as the non-polar analog of similar size and shape as the species of interest. Unfortunately this method breaks down when it comes to ringed structures and halogenated species. Quite often the polar parameter is calculated from Hansen and Beerbower’s simplification of the Böttcher equation (12): 37.4 / / (5) Where  is the dipole moment in Debye and VM is the molar volume in cm3/mol. This requires the user to know or estimate these values with some degree of certainty. When this cannot be done, group contribution methods are usually employed. When the dispersion and polar solubility parameters can be estimated with confidence the hydrogen-bonding parameter can be calculated from equation 4; this requires the use of energy of vaporization data to determine the total solubility parameter. When either P or D cannot be estimated to solve for H, then the group contribution method outlined by Stefanis and Panayioutou appear to be sufficient to calculate the hydrogen-bonding parameter independently. 8 Current Project The funding agency for this project, B&W Pantex in Amarillo, TX is interested in the solubility of nitrated aromatic and nitramine species such as 1,3,5-Trinitro-2-[2(2,4,6-trinitrophenyl)ethenyl]benzene (HNS) and 1,3,5,7-tetranitro-1,3,5,7-tetrazocane (HMX), Figure 1. HNS is an example of a nitrated aromatic species and HMX is an example of a nitramine. These species exhibit poor solubility in common industrial solvents. The best solvents for many of these applications are N-methylpyrrolidone and N,N-dimethylformamide though they do not perform as well as desired. The Hansen solubility parameters of these explosives and other energetic materials will be explored to determine if these two solvents are the best solvents, or if blended solvents are better for dissolving energetic materials. NO2 NO2 O 2N O 2N N I. NO2 O 2N N II. N NO2 N NO2 O 2N NO2 Figure 1. Chemical Structures of HNS (I) and HMX (II). There are significant problems with the estimation of the HSPs of nitrated aromatics and nitramines. Primarily, they have poor solubility in most common industrial solvents. This precludes a wet chemistry solubility study. Secondly, there are no entries in group contribution methods for aromatic nitro-groups or nitramines. This precludes all group contribution methods for the estimation of HSPs. The Yamamoto 9 group contribution method present in HSPiP contains nitro groups, but gives unreasonable results, and is still under development. One could approach this problem with a brute force technique and try various combinations of two, three, and four component solvent blends to determine which matches the HSPs of the unknown. The number of combinations of two, three, and four component blends of a possible ten solvents is 375. If each combination is tested at three different volume fractions, the number of experiments top one thousand. If the number of possible solvents is expanded by just one more to eleven, then the number of experiments would top two thousand. In this case, time and cost preclude the use of this method for the determination of HSPs. Problems of this size and complexity often necessitate the use of computers. It is the aim of this research to simplify and accelerate the determination of HSPs through the use of electronic structure calculations. These calculations are cheap, relatively accurate, and if employed with modern parallel processing computers they can be fast. Structure variables will be extracted from calculation result files and used to develop a quantitative structure activity relationship (QSAR). Electronic structure calculations will provide the structure variables and a matrix of coefficients will then relate them to HSPs (the desired activity). The result will be an inexpensive and unambiguous method for calculating the HSPs of nitrated aromatics and nitramines. Theory As recent as the 1960s, chemists have begun to use computers to aid them in the modeling of molecular structure (13). There have been many developments since then and many programs now are capable of the accurate modeling of bond lengths, molecular 10 orbitals, and thermodynamic properties. These computational chemistry models allow scientists to predict the bulk properties of a species from the calculation of a single molecule. There are significant advantages to this. One advantage is cost. If a scientist can research a species or a group of species without ever having to synthesize them, then the cost of conducting that research is tremendously reduced. Aside from the obvious cost of consumable materials required to synthesize a species there are many less tangible costs. For example, the user’s time spent on synthesis, purification, characterization, and evaluation. Reduced instrument maintenance costs as well as less exposure to laboratory chemicals are also seen as benefits. The least computationally expensive (fastest) model chemistry is molecular mechanics which uses pair potentials between all atoms in the system. These potentials have been parameterized to match the most common types of molecules. Unfortunately the highly nitrated species of interest are not well targeted by the molecular mechanics parameter set. Electronic structure calculations seek to account for the electrons’ contribution to the chemistry of a species. Because the exact location and velocity of an electron cannot be known, the numerical Schrödinger equation cannot be solved explicitly. The pieces which cannot be solved must be approximated by various levels of electron-correlation theory. Most structure calculations require a designation of level of theory and basis set of wavefunctions. Some popular levels of theory include HF (Hartree-Fock) and B3LYP (Becke type, 3-parameter exchange with Lee-Yang-Parr correlation functional) (14). The basis set defines the number and type of functions that will be used to approximate the 11 wavefunction of each electron in the species. A typical basis set is 6-31G(d). This nomenclature indicates that six functions will be used to approximate the wavefunction of the core electrons, three functions will be used for the valence electrons of heavy atoms (all atoms which are not hydrogen), and one function will be used for hydrogens. The “G” indicates that the type of functions to be used will be Gaussian functions; and the “(d)” indicates that d-orbital polarization functions will be added to heavy atoms. Typically the input geometry is not the lowest energy conformation and so the molecule will require optimization. Any number of subsequent calculations can be performed on the resulting optimized geometry. The most common calculation is a “frequency calculation.” Many of the electronic structure parameters are calculated in this step, the most prominent of which are the vibrational frequencies, and infrared(IR) and Raman spectral intensities. The publication syntax for the specification of a geometry optimization with Hartree-Fock theory and a 3-21G(d) basis set followed by a frequency calculation of higher level of theory and basis set would be: B3LYP/631G(d,p)//HF/3-21G(d). The B3LYP/6-31G(d) model chemistry was chosen for this project because of its high accuracy and its low computational cost. The published mean average deviation for this method on a standard set of molecules is 7.9 kcal/mol with a standard deviation of 9.5 kcal/mol (14). The utility of electronic structure calculations to approximate bulk properties via a structure activity relationship has been demonstrated (15). This technique employs the use of structure variables to predict an unknown activity. Each variable is assigned a coefficient which indicates its weight or effect toward the total activity. For a quantitative structure activity relationship (QSAR) to be established the relationship 12 matrix must be trained to a set of species for which the activity and structure variables are well defined. This will yield a matrix which can be applied to the prediction of activity for species which are sufficiently similar to the training set. This relationship “training” can be accomplished by a step-wise regression technique. The step-wise regression is a methodical way for including or excluding structure variables from the relationship matrix. The driving force for the inclusion or exclusion of variables is their correlation with activity. Correlation is quantified in the form of an Rsquared value. The criteria for the decision to include or exclude a particular variable is the probability that random noise could produce the coefficient of that variable in the structure-activity transformation matrix. This criteria is quantified by a P-test which results in the appropriately named P-value. A variable is chosen and a regression is performed, producing a transformation coefficient. There will be some level of correlation with this first structure variable. Another variable is included and another regression is performed. The coefficients of the two variables and their P-values are compared to the exclusion criteria. If a third variable is inserted and causes the P-value of another variable coefficient already in the model to exceed the exclusion criteria and also causes an increase in the R-squared value, then the third variable is kept and the criteria-exceeding variable is ejected from the model. This process is repeated until all variables in the model meet the inclusion criteria. In the case of this research each Hansen solubility parameter will have a relationship matrix and thus will require a step-wise regression. This means that for every model there will be three step-wise regressions on n number of variables which will yield three 1 × n matrices ([D]n, [P]n, and [H]n) or one 3 × n matrix ([D P H]n). These 13 matrices can be simplified by eliminating those variables which do not have a statistical contribution to the activity (coefficient = 0) to yield matrices with numbers of rows less than n.