A general purpose acetonitrile interaction potential to describe its liquid, solid and gas phases

Graphical Abstract1 A General Purpose Acetonitrile Interaction Potential to Describe its Liquid, Solid and Gas2 Phases3 Jorge Hernández-Cobos,José M. Martínez,Rafael R. Pappalardo,Iván Ortega-Blake,Enrique Sánchez Marcos4 5 Highlights1 A General Purpose Acetonitrile Interaction Potential to Describe its Liquid, Solid and Gas2 Phases3 Jorge Hernández-Cobos,José M. Martínez,Rafael R. Pappalardo,Iván Ortega-Blake,Enrique Sánchez Marcos4 •Acetonitrile flexible and polarizable force field based on quantum mechanical computations5 •MD simulations of gas phase clusters, liquid, and 𝛼and 𝛽crystal structures6 •Satisfactory comparison of simulated and experimental physicochemical properties7 A General Purpose Acetonitrile Interaction Potential to Describe its1 Liquid, Solid and Gas Phases⋆ 2 Jorge Hernández-Cobosa,José M. Martínezb,Rafael R. Pappalardob,Iván Ortega-Blakeaand3 Enrique Sánchez Marcosb 4 aInstituto de Ciencias Físicas, UNAM, Apartado Postal 48-3, 62251-Cuernavaca, México5 bDepartament of Physical Chemistry, University of Sevilla, 41012-Sevilla, Spain6 7 ARTICLE INFO Keywords: ab initio force field polarizable and flexible acetonitrile model MD simulations x-ray and neutron diffractograms physicochemical properties acetonitrile crystals 8ABSTRACT9 10 A flexible and polarizable force field to describe acetonitrile in its three states of matter has11 been developed on the basis of an ab initio potential energy surface at the MP2 level. The12 acetonitrile molecule is represented by a sophisticated twelve-site model: the six nuclei plus13 six mobile charges. Three thousands structures have been employed for the fitting including14 monomers, dimers and trimers. Gas phase behavior has been tested by analyzing aggregates15 from dimers to 27-mer which give accurate reproducibility of the experimental and quantum-16 mechanical results. Classical MD simulations have been performed for the liquid and its17 crystalline 𝛼−and 𝛽−forms. An analysis of liquid in a wide range of structural, dynamic,18 energetics and spectroscopic properties has been performed, as well as neutron diffraction data.19 Comparison with experimental information is satisfactory. Due to the first principles nature and20 the polarizable character of this force field properties of the three states of matter are properly21 reproduced.22 23 1. Introduction24 Acetonitrile (ACN) is one of the most used polar aprotic solvents, more precisely a dipolar non-hydrogen bonded25 solvent.[1] ACN combines a high relative permittivity, low viscosity, differential solvation behavior against cations26 and anions, and an amphiphillic character provided by the presence of the methyl group and a strong dipolar charge27 distribution of the molecule. The concurrence of all these features makes this solvent extremely versatile in manifold28 technological areas such as chemical and pharmaceutical industries,[2,3], electrochemistry and organic and inorganic29 synthesis.[1,4] Even more, ACN plays a double role as solvent and reactant in solution as well as in the acetonitrilation30 of gaseous aggregates.[5,6]31 A thorough understanding of the physicochemical properties of liquid ACN and its solutions needs their structural32 characterization, not only to explain the observed behavior but also to envisage the study of more complex chemical33 systems,andthe extensionto newapplications ingasphase[6] andsolid state.An excellentreview ofPothoczkietal.[7],34 appeared in 2015, compiles the experimental and theoretical efforts to deal with liquid ACN. The set of experimental35 x-ray and neutron diffraction studies[8–12] has been complemented in Table 1with several theoretical methods to get36 insight into the structural interpretation at the atomistic level.[13–25]37 The number of particles needed to describe the condensed medium, the performance of the available DFT38 methods, as well as the long simulation times required to estimate certain properties limit the use of ab initio MD39 simulations.[26,27] Then, most of the statistical simulations are classical and the treatment of the condensed phase40 relies on the force fields describing the interactions among particles. Different strategies have been adopted to provide41 good-behaving interaction potentials, considering on one hand, the flexibility or rigidity of the ACN molecule and42 on the other the polarizable character of the molecule. Likewise, the data used to build the force field could be43 from experimental or first principles origin. In the first case, the potential parameters are fitted to reproduce as44 well as possible some selected experimental properties. In the second case, the fitting is based on the reproduction45 of quantum mechanical results of the ACN monomer and, in general, the dimer interaction in a large number of46 attractive and repulsive arrangements. Table 1collects the force fields available in the literature with indication of47 their main features and computed properties. Most of them are able to supply good descriptions of a set of ACN liquid48 ∗Corresponding author [email protected] (E. Sánchez Marcos) ORCID(s): 0000-0002-8367-9105 (E. Sánchez Marcos) :Accepted manuscript Page 1 of 18 Simulation of acetonitrile in liquid, gas and solid phases properties whose range of applicability is determined by both, the type of model defined and the set of data used for1 the fitting. In a recent and clarifying study, Kowsari et al.[13] have performed a systematic evaluation of the ACN2 force fields which use a six-site model and molecular flexibility.[14,15,17,20] They compared on the same foot the3 performance of these potentials on a set of thermodynamic and transport properties of liquid ACN. These authors point4 out the difficulty of satisfying simultaneously very accurate description of thermodynamic and transport properties.5 Thus, the Grabuleda et al.’s potential[20] and that of Nikitin et al.[17] supply good description of the density of the6 vaporization enthalpy, but they provide less accurate values of the self-diffusion coefficient. Contrarily, Orhan[15]7 builds a force field where the transport properties, the dielectric relaxation and vibrational spectra of bulk ACN were8 wellreproduced along awide temperaturerange.However, Kowsarietal.[13] showthatthedensity andthe vaporization9 enthalpy are underestimated by Orhan’s potential. Koverga et al.[14] adopt an intermediate strategy to develop their10 force field since they fit parameters to reproduce simultaneously thermodynamic and transport properties. They get a11 good potential which is able to predict properly density, vaporization enthalpy and diffusion coefficient compared to12 the rest of available simulated data. Nevertheless, overall the dynamic properties show a lower accuracy relative to13 thermodynamic estimates.[14] On the basis of this analysis, Kowsari et al.[13] performed a refinement process tuning14 some of the Lennard-Jones parameters of the potentials to fit the available experimental data. The refined potential is15 able to reproduce both thermodynamic and dynamic properties of liquid ACN pretty well.16 It must be noted, that there is no purely first-principles interaction potential including polarizability and flexibility17 on a fully atomic representation of the ACN molecule. In this work we do build such a force field type, in addition18 of including many-body corrections. To this end we have developed a twelve-site model.The aim is to check the19 performance of this purely first-principles force field, not only in liquid state but also in gas phase clusters and solid20 state crystals. Since no experimental information is used for the force field building, the inclusion of enough many-21 body contributions should allow a proper description of all the aggregation states involving ACN molecules. The gas22 phase was studied on ACN clusters of different sizes and the solid state in two well-characterized crystal forms, 𝛼and23 𝛽.[28–30]24 Table 1: Computer simulation studies of ACN based on classical interaction potentials. Author[Ref.] Potential Nature Geometry Charge Distrib. Potential Fitted to(𝑎) Properties(𝑎) Kowsari [13] Flexible Non-Polarizable Exp. data of 𝜌, 𝑉m,Δ𝐻vap, , C.E.D. and D 𝜌,𝑉m,Δ𝐻vap, C.E.D., D, RDF, CDF Koverga [14] Flexible Non-Polarizable CPMD, QM, Exp. data RDF, CDF, Nearest neighbors analysis, Dielectric properties, H-bond interactions Orhan [15] Flexible Non-Polarizable Exp. data RDF, Liquid vibration spectra at different T, 𝜂s,𝐷,𝜏1, 𝜏2, Dielectric properties Alberti [16] Rigid Non-polarizable Atom polarizabilities Econfig as a function of T (CH3CN)n, n=2-4 𝐷and 𝜌as a function of T RDF, Ion-ACN solutions Nikitin [17] Flexible Non Polarizable Exp. data of 𝜌, Δ𝐻vap RDF, 𝜌of ACN-water solutions, Dielectric properties, 𝐷 Gee [18] Flexible Non-polarizable Exp. data Δ𝐻vap,𝜌, CP,𝜀,𝐷,𝜂,𝜏1, 𝜏2,RDF :Accepted manuscript Page 2 of 18 Simulation of acetonitrile in liquid, gas and solid phases Author[Ref.]/Pot. Name Potential Nature Geometry Charge Distrib. Fitted to Properties Sum [31] Rigid Non-Polarizable & Polarizable Scaled ab initio 2nd. Virial Coefficient Δ𝐻vap, Gas-Liquid Coexistence: Tc,𝜌c RDF Guardia [19] Rigid Non-Polarizable ab initio calculations of dimer 𝜌,𝑉m,ΔHvap,𝜅, CpRDF, X- ray structure factor, 𝐷,𝜏1, 𝜏2Dielectric properties (Gk, y, 𝜀0,𝜀o-𝜀∞, gk,𝜏𝜑,𝜏D, 𝜏D/𝜏) Grabuleda [20] Flexible Non-polarizable Exp. data 𝜌,𝜅,Δ𝐻vap 𝜇, RDF Cabaleiro-Lago [21] Rigid Polarizable ab initio IMPT data HF/MP2/ 6-31G** Dimer energy, geometry and frequencies Trimer and tetramer energy and geometry 𝑈, RDF, Dipole moment distribution, Dipole-Dipole angle distribution Siebers [22] Flexible Non-Polarizable Exp. data Minimum energy structures of ACN clusters from n = 2 to n = 13 Frequency shifts of the CC stretch and the CH3wag modes Jorgensen [23] Rigid Non-Polarizable Exp. Data 𝜌,𝑉M,𝑈,ΔHvap,𝜅, Cp, RDF, Energy distribution functions, RDF, Dipole correlation Edwards [24] Rigid Non-Polarizable Six site model Dielectric Properties Center of mass G(r) Orientational correlation functions Dipole correlation functions Böhm [25] Rigid Non-Polarizable Exp. data 𝑈,𝜅, Cp,RDF, 𝐷,𝜏2Neu- tron and X-Ray Scattering (𝑎)𝜌= density; 𝑉m= molar volume; C.E.D.= Cohesive energy density;1 CDF=Combined distribution function; 𝜂= shear viscosity; 𝐷= self-diffusion coefficient;2 𝜏i= i𝑡ℎ-order reorientational time; 𝐶p=heat capacity at constant pressure; 𝜅= isothermal compressibility.3 2. Methods4 2.1. Potential Development5 In this work we build a force field for acetonitrile whose functional form is based on the Mobile Charge Densities6 in Harmonic Oscillators (MCDHO) model [32,33]. This allows for the inclusion of intramolecular flexibility, non-7 additivity and polarizability. The molecular polarizability is accounted by a set of negative mobile charge densities8 with radial exponential decay associated to each positively charged nuclei. The interaction between both charges is9 :Accepted manuscript Page 3 of 18 Simulation of acetonitrile in liquid, gas and solid phases Figure 1: Molecular model used to define the acetonitrile molecule. described by means of a harmonic oscillator. The remaining charge-charge interactions in the system consider the1 charges as points. A mobile charge has been associated to each ACN nucleus, thus we build a twelve-site model formed2 by the sum of six nuclei plus six mobile charges. Figure 1displays a scheme of the ACN molecular model.3 The specific expression for the water - water interaction in the MCDHO model was given in eqs. (4-6) of Ref.[32].4 The molecular complexity of acetonitrile requires a more general formulation that we detail at this point. A previous5 force field, developed for liquid methanol [34], uses a model definition similar to that here presented. However, in that6 case hard core cutoff radii for the interactions were included in the potential limiting its use to Monte Carlo simulations.7 For the ACN we developed a potential with a good behavior of the energy gradients, allowing for its use in Molecular8 Dynamics simulations.9 The total energy is the sum of an intramolecular and an intermolecular energy. In the following equations10 𝑟𝑖𝑗 =|𝐫𝑖−𝐫𝑗|is the distance between the corresponding centers, either fixed positive nuclei 𝑍or mobile negative11 charges 𝑞, where 𝑖and 𝑗subscripts correspond to sites 𝑖and 𝑗.12 The analytical expression to compute the intramolecular energy is the sum of the following contributions:13 𝑈intra =∑ 𝑖∈𝑆∑ 𝑗>𝑖 [𝑈(𝑍𝑖, 𝑍𝑗) + 𝑈(𝑞𝑖, 𝑞𝑗)]+∑ 𝑖∈𝑆∑ 𝑗≠𝑖 𝑈(𝑞𝑖, 𝑍𝑗) + ∑ 𝑖∈𝑆 1 2𝑘𝑖𝑟2 𝑖𝑖 (1) +∑ bonds 𝑈𝑏(𝑖, 𝑗) + ∑ angles 𝑈𝜃(𝑖, 𝑗, 𝑘) where the first two terms correspond to the electrostatic interactions among nuclei and mobile charges, the third term14 corresponds to the harmonic energy of the mobile charge at a distance 𝑟𝑖𝑖 from its nucleus, and the other two terms15 correspond to the bond energy and the bond angle energy. The summation 𝑖∈𝑆runs over atomic centers of the16 molecule. These terms do have the following functional forms:17 1. The electrostatic interaction between nuclei:18 𝑈(𝑍𝑖, 𝑍𝑗) = 𝑍𝑖𝑍𝑗 𝑟𝑖𝑗 (2) 2. The interaction between mobile charges with decay lengths 𝜆𝑖and 𝜆𝑗attached to bonded atoms. In this case a19 two-center integral should be used; however, an approximate expression that gives good accuracy at all relevant20 distances was found[32,33]:21 𝑈(𝑞𝑖, 𝑞𝑗) = 𝑞𝑖𝑞𝑗 𝑟𝑖𝑗 ⎡⎢⎢⎣ 1 − ⎛⎜⎜⎝ 𝜆𝑖𝜆𝑗+(𝜆𝑖+𝜆𝑗)2 (𝜆𝑖+𝜆𝑗)3𝑟𝑖𝑗 + 1⎞⎟⎟⎠ exp ⎛⎜⎜⎝ −2𝜆𝑖𝜆𝑗+(𝜆𝑖+𝜆𝑗)2 (𝜆𝑖+𝜆𝑗)3𝑟𝑖𝑗⎞⎟⎟⎠⎤⎥⎥⎦ (3) mobile charge densities are considered as point charges when interacting with mobile charge densities of non-22 bonded atoms23 𝑈(𝑞𝑖, 𝑞𝑗) = 𝑞𝑖𝑞𝑗 𝑟𝑖𝑗 (4) :Accepted manuscript Page 4 of 18 Simulation of acetonitrile in liquid, gas and solid phases 3. The electrostatic interaction between mobile charge densities 𝑞𝑖and an atomic center charge 𝑍𝑗to which its own1 atom is bonded, where 𝜆𝑖is the intramolecular decay length of the charge density:2 𝑈(𝑞𝑖, 𝑍𝑗) = 𝑞𝑖𝑍𝑗 𝑟𝑖𝑗 [1 − (𝑟𝑖𝑗 𝜆𝑖 + 1)exp (−2𝑟𝑖𝑗 𝜆𝑖)] (5) The mobile charge densities are considered as point charges when interacting with a non-bonded atomic center.3 In these cases the electrostatic interaction is given by:4 𝑈(𝑞𝑖, 𝑍𝑗) = 𝑞𝑖𝑍𝑗 𝑟𝑖𝑗 (6) 4. A Morse potential between each pair of bonded atoms, with depth 𝐷𝑖𝑗 , inverse decay length 𝛾𝑖𝑗 and equilibrium5 parameter 𝑟𝑒𝑞 𝑖𝑗 :6 𝑈𝑏(𝑖, 𝑗) = 𝐷𝑖𝑗 {exp [−2𝛾𝑖𝑗 (𝑟𝑖𝑗 −𝑟𝑒𝑞 𝑖𝑗 )]− exp [−𝛾𝑖𝑗 (𝑟𝑖𝑗 −𝑟𝑒𝑞 𝑖𝑗 )]} (7) 5. A quadratic term for each bond angle 𝜃𝑖𝑗𝑘, with parameter 𝑘𝜃and the equilibrium value 𝜃0 𝑖𝑗𝑘:7 𝑈𝜃(𝑖, 𝑗, 𝑘) = 𝑘𝜃 2(𝜃𝑖𝑗𝑘 −𝜃0 𝑖𝑗𝑘)2(8) The intermolecular energy is composed of the following terms:8 𝑈inter = 𝑁 ∑ 𝑆=1 (𝑆−1 ∑ 𝑇=1 {∑ 𝑖∈𝑆∑ 𝑗∈𝑇[𝑈inter(𝑍𝑖, 𝑍𝑗) +𝑈inter(𝑞𝑖, 𝑍𝑗) + 𝑈inter(𝑞𝑗, 𝑍𝑖) + 𝑈inter(𝑞𝑖, 𝑞𝑗)]}) (9) where summations 𝑖∈𝑆and 𝑗∈𝑇refer to entire molecules and the expression of each term is:9 1. The electrostatic interactions between atomic centers plus exponential terms with parameters 𝐴𝑖𝑗,𝐵𝑖𝑗,𝑎𝑖𝑗 and10 𝑏𝑖𝑗;11 𝑈inter(𝑍𝑖, 𝑍𝑗) = 𝐴𝑖𝑗 𝑒𝑥𝑝(−𝑎𝑖𝑗 𝑟𝑖𝑗 ) + 𝐵𝑖𝑗 𝑒𝑥𝑝(−𝑏𝑖𝑗 𝑟𝑖𝑗 ) + 𝑍𝑖𝑍𝑗 𝑟𝑖𝑗 (10) 2. The electrostatic interaction between mobile charges 𝑞𝑖considered as point charges and charges 𝑍𝑗:12 𝑈inter(𝑞𝑖, 𝑍𝑗) = 𝑞𝑖𝑍𝑗 𝑟𝑖𝑗 (11) 3. The electrostatic interaction between mobile charges,13 𝑈inter(𝑞𝑖, 𝑞𝑗) = 𝑞𝑖𝑞𝑗 𝑟𝑖𝑗 ,(12) :Accepted manuscript Page 5 of 18 Simulation of acetonitrile in liquid, gas and solid phases Therefore, the energy of a cluster with 𝑁molecules is given by,1 𝑈total =𝑈inter + 𝑁 ∑ 𝑆=1 𝑈intra(𝑆),(13) The interaction energy, 𝑈int, for a cluster of N molecules requires the subtraction of the intramolecular energies of2 the isolated molecules, 𝑈0 intra, that corresponds to the 𝑈intra value in eq. (1) for an ACN molecule at its optimized3 geometry.4 𝑈int =𝑈total − 𝑁 ∑ 𝑆=1 𝑈0 intra(𝑆) = 𝑈total −𝑁𝑈0 intra,(14) thus, taking into consideration the energetic cost of polarizing and deforming each molecule in the cluster or condensed5 phases.6 Due to the large number of parameters to be fitted in this potential and the different nature of the interactions,7 the set of ACN structures chosen to be fitted from the ab initio potential energy surfaces were chosen in an iterative8 manner. Likewise, the different set of parameters were included sequentially in the fitting process as well. That is,9 initial structures were taken from a regular scan of the ACN molecule hyper-surface with respect to their different10 internal degrees of freedom. An initial potential parameterization was fitted for the intramolecular part of the potential,11 𝑈intra. This potential was used to generate a set of new structures derived from numerical simulations at standard12 temperature. These new structures were then computed at the ab initio level and added to the set of structures employed13 to fit the potential. This procedure is repeated until convergence is achieved. A scheme of this iterative process was14 sketched in a previous work where a Cf(III)-H2O potential was built.[35]. Once a fitting of the intramolecular potential15 was achieved, ACN dimers were included in the energy surface, scanning intermolecular distances and orientations.16 Then, the same iterative procedure was used for the intermolecular part adding dimer interactions energies and non-17 additivities computed from trimers. Hence, the final surfaces to be fitted consisted of: (a) the three components of18 the dipole moment of the optimal monomer and the dipole moment modulus of all distorted ACN monomers. (b)19 1117 monomer structures where the deformation energy, relative to the optimal monomer, was considered for fitting,20 (c) 1980 dimer structures where the interaction energy was fitted, (d) 61 trimer structures where the three-body non21 additive contributions were fitted.22 The quality of the fitting is shown in Figure 2where the total interaction energy, 𝑈int in eq. (14), computed23 by the force field for the employed set of 3158 structures is compared with the corresponding ab initio values. The24 standard deviation of this comparison is 0.29 kcal/mol. Due to the fitting procedure four different standard deviation25 can be evaluated: that of the monomer deformation energies, 0.25 kcal/mol, that of the dimer interaction energies, 0.2826 kcal/mol, that of the many-body interactions associated to the trimers, 0.04 kcal/mol, and that of the dipole moment27 evaluation of monomers, 0.11 D. The potential parameters are presented in Table S1 of Supplementary Material (SM).28 Molecular orbital calculations were done with the Gaussian09 program [36]. The computational level chosen29 was MP2/6-311++G(3df,3pd). Interaction energies were corrected of basis set superposition errors by using the30 counterpoise method.31 Table 2collects the geometrical parameters corresponding to an ACN molecule and its most stable dimer32 obtained by the quantum-mechanical calculation, the developed potential, results obtained by other authors and33 experimental data. Both the degree of reproducibility of the basic molecular properties and intermolecular interactions34 are satisfactory. Table S2 of SM collects the comparison of the frequencies of the main vibrational normal modes of35 the ACN molecule obtained by the force field with those derived from the QM level used for the fitting and with the36 experimental ones. More demanding tests to assess the good performance of the developed potential, when a large37 number of ACN molecules are interacting, are presented in the Results Section.38 2.2. Molecular Dynamics Simulation Details39 Molecular dynamics (MD) simulations were run using a modified version of DL_POLY_CLASSIC [45] which40 allows the use of the MCDHO force field. An integration time step of 0.05 fs was used to account properly for the41 high frequencies associated to the dynamics of the mobile charges. The Ewald Summation was used to calculate the42 :Accepted manuscript Page 6 of 18 Simulation of acetonitrile in liquid, gas and solid phases Figure 2: Comparison between the interaction energies, 𝑈int predicted by the force field developed and the ab initio ones for the set of structures employed for the fit. Table 2 Basic properties of ACN molecule and its dimer. Property QM Potential Literature𝑎Exptal.𝑏 monomer 𝑅C-C (Å) 1.458 1.457 1.441-1.464 1.458 𝑅C-N (Å) 1.168 1.170 1.159-1.181 1.157 𝑅C-H (Å) 1.087 1.087 1.093-1.109 1.104 ΘCCH (◦) 109.0 109.0 109.0-110.5 109.5 𝜇(D) 4.3 4.0 3.98-4.12 3.9 dimer 𝑅(CN...CN) (Å) 3.291 3.369 3.096-3.502 𝑅(N...H) (Å) 2.468 2.563 2.450-2.601 Θ(N-C...CN) (Å) 134.9 134.3 132.9-139.3 𝑈int (kcal/mol) -6.21 -6.22 (-3.08)-(-9.85) -5.0 𝑎Refs. [16,37–39]𝑏Refs. [40–44] electrostatic interactions. The van der Waals interactions were truncated at a cutoff radius of 14 Å. The Nosé-Hoover1 thermostat to attain the average target temperature was applied with a characteristic time of 0.5 ps. The Nosé-Hoover2 barostat with a characteristic time of 1.0 ps was applied to keep the average pressure at 1 atm. 500 ACN molecules were3 explicitly considered in the simulation box for the liquid ACN simulation at 298 K. The system was minimized and an4 NPT run of 0.5 ns to equilibrate the simulation box was performed. 1 ns production simulation was used for further5 analysis and extraction of physicochemical properties. Additionally, a similar NVT MD simulation at the experimental6 density at 298 K was also performed to compare with the NPT results. The translational self-diffusion coefficient, 𝐷0,7 for ACN in the MD simulation was computed using the Einstein formula[46] and corrected for finite box-size effects8 with the Yeh and Hummer’formula [47]:9 𝐷corr ACN =𝐷0 ACN +𝑘𝐵𝑇 𝜉 6𝜋𝜂𝐿 (15) :Accepted manuscript Page 7 of 18 Simulation of acetonitrile in liquid, gas and solid phases Table 3 Basic properties of liquid ACN at 298 K. Property𝑜MD simulation Literature Exptal. Δ𝐻vap (kJ mol−1) -33.0 (-30.7)-(-34.2)𝑎(-32.9)-(-34.5)𝑏 𝜌(kg m−3) 793 718-776𝑐776-778𝑑 𝜅(GPa−1) 1.07 1.18-1.62𝑒1.07-1.13𝑓 𝜇(D) 4.67 4.01-4.88𝑔4.5ℎ 𝐷·109(m2s−1) 3.81 3.02-4.93𝑖4.04-5.10𝑗 𝜏1(ps) 3.9 3.28-3.68𝑘 𝜏2(ps) 1.3 1.02𝑙 𝜂(cP) 0.43 0.29-0.64𝑚0.34𝑛 𝑎Refs.[14,15,17,20,52,53]𝑏Refs.[54,55]𝑐Refs.[14,15,17,20,27,52,53] 𝑑Refs.[13,55,56]𝑒Refs.[52]𝑓Refs.[57]𝑔Refs.[20,21,58]ℎRefs.[59] 𝑖Refs.[14,15,17,20,52,53]𝑗compilation in Ref.[13]𝑘Refs.[60,61]𝑚Ref.[18,19,23,24] 𝑛Ref.[62] 𝑜vaporization enthalpy, Δ𝐻vap; density, 𝜌; isothermal compressibility, 𝜅; dipole moment, 𝜇; diffusion coefficient, 𝐷; 1st. and 2nd. dipole moment axis reorientational times, 𝜏1and 𝜏2;𝜂, viscosity. Figure 11: Neutron diffraction data for liquid ACN at 293K: experiment (black dots) taken from ref. [8], simulated data from NPT MD simulation(red line) and from NVT MD simulation with the experimental density (blue line) the slight density change between the NPT and NVT simulated functions. It is worthy to comment that Pothoczki and1 Pusztai[12] using the six-site OPLS force field of Jorgenssen and col.[64] computed the X-ray and neutron diffraction2 data of ACN obtaining an excellent agreement with the experimental data.3 Since the developed force field uses a flexible and polarizable molecular model for ACN, a final analysis has been4 the inspection of how the average dipole moment of ACN molecules, < 𝜇 >, changes for the different phases. This5 evolution is presented in Figure 12 where the mean dipole moment for the set of gas phase clusters previously examined6 in Figure 5is plotted as a function of the number of molecules forming the aggregate. The vertical bars correspond to7 the dispersion of the individual dipole moments present in the clusters considered. It should be noted that the < 𝜇 >8 increases with the cluster converging to a value close to that of the liquid MD simulation. This is not the case for the9 27-mer cluster that presents an average dipole moment value larger than the mean liquid value. This is due to the fact10 :Accepted manuscript Page 14 of 18 Simulation of acetonitrile in liquid, gas and solid phases Figure 12: Evolution of the average dipole moment (D) of ACN molecule with the cluster size, indigo dot corresponds to a (ACN)27 cluster with the structure corresponding to the ACN crystal arrangement. Vertical bars indicates the range of dipole values of molecules in a given cluster. Horizontal dashed lines indicates the mean value of dipole moment in liquid at 300K (red) and in the 𝛼(blue) and 𝛽crystals (green). that the structure of this cluster is a fragment of the 𝛼-form crystal structure, where dipole alignment is predominant1 (see Figure 6). This fact is confirmed by the computed mean value of the crystal 𝛼-form (blue line) which is larger than2 that of the liquid and the 27-mer cluster. Finally the structural arrangement of the 𝛽-form leads to a higher average3 dipole moment. Apart from an intrinsic good description of the charge distribution provided by the force field for the4 isolated molecule yielding 4.0 D vs. an experimental value of 3.96 D,[42] the molecular polarization of our model is5 able to reproduce environmental effects progressively more complex by including many-body and temperature effects.6 4. Concluding remarks7 This study has shown the validity of a new flexible and polarizable force field to deal with the three states of8 matter of acetonitrile based on first principles calculations exclusively. It has been shown that the MCDHO strategy9 is versatile enough to grasp the polarizable character of acetonitrile under different environments, going from small10 molecular aggregates to either liquid at room temperature or two crystaline forms. Hence this study highlights the11 importance of a good description of both, the electronic polarization and the intramolecular flexibility together with12 the explicit inclusion of many-body contributions. This permits the proper reproduction of the molecular and bulk13 behavior at different phases. It is also shown that a proper molecular description of the acetonitrile can be reached14 without including experimental information in the force field building, although a sophisticated formalism was needed15 to develop it.16 The consistency of the potential supports two conclusions. Contrary to the case of the ubiquous solvent, water,17 quantum-mechanical ingredients incorporated through the MP2method with alarge even-temperedbasis setsis enough18 to describe both the molecular and the condensed phase behavior of acetonitrile. The second conclusion is that many-19 body effects can be incorporated by explicit molecular polarizability and a limited number of acetonitrile molecules20 in the clusters considered in the ab initio potential surfaces used for the fitting.21 The structural analysis of liquid ACN reveals that its structure extends to the first solvation shell where a dominant22 pair-orientation does not exits, antiparallel, cross, H-bond and T1-shaped being the more frequent arrangements as a23 function of the intermolecular distance.24 :Accepted manuscript Page 15 of 18 Simulation of acetonitrile in liquid, gas and solid phases The classical nature of the potential will allow the exploration of highly-demanding computational properties of1 acetonitrile liquid, as dielectric and rheological properties, because of the large number of molecules to be considered2 or the need of long simulation times are affordable. The high-level ab initio nature of the information employed to3 characterize the ACN molecule and its interactions wich other molecules under many different arrangements allows a4 refined and unbiassed molecular answer in different environments.5 The combination of the two previous feature potential leads to the conclusion that this potential gives an even-6 tempereddescription ofACN-containingsystemswherea compromisebetweencondensed phasefactors andmolecular7 properties simultaneously contribute to the physicochemical porperties of the system. Then future efforts will be8 addressed to find out the limits of the force field by computing significant bulk-depending mechanical, dynamical,9 electric, magnetic , rheological and spectroscopic properties of acetonitrile under a wide range of temperature and10 pressure. In addition, this potential could be applied in multicomponent systems in different agreggation states: gas,11 solutions, confined media, interphases and solid and crystalline systems.12 5. Acknowledgements13 We thank Dr. Silvia Pothoczki for providing us the experimental neutron diffraction data of liquid acetonitrile.14 We are grateful to Reyes Garcia Carreon for support in the computational part of this work. This document is the15 result of the research project (PGC2018-099366-B-I00) funded by the Spanish Ministry of Science, Innovation and16 Universities and by the Universidad Nacional Autónoma de México(UNAM)(DGAPA-PAPIIT-IG100920) and Project17 LANCAD-UNAM-DGTIC-057 for computational resources.18 6. Appendix A. Supplementary data19 Supplementary data to this article can be found online at20 References21 [1] C. Reichardt, Solvent and Solvent Effects in Organic Chemistry, 2nd ed., Wiley-VCH, Weinheim, 2004.22 [2] N. Kumar, D.Sangeetha, P. Balakrishna, Development and validation of a UPLC method for the determination of duloxetine hydrochloride23 residues on pharmaceutical manufacturing equipment surfaces, Pharm. Methods 2 (2010) 161–166.24 [3] L. Wannatong, A. Sirivat, P. Supaphol, Effects of solvents on electrospun polymeric fibers: preliminary study on polystyrene, Pharm. Methods25 53 (2004) 1851–1859.26 [4] J. Kim, D.-J. Lee, H.-G. Jung, Y.-K. Sun, J. Hassoun, B. Scrosati, Lithium-sulfur batteries: an advanced lithium-sulfur battery, Adv. Funct.27 Matter. 23 (2013) 1076–1080.28 [5] G. Caldwell, M. D. Rozeboom, J. P. Kiplinger, J. E. Bartmess, Displacement, proton transfer, or hydrolysis? mechanistic control of acetonitrile29 reactivity by stepwise solvation of reactants, J. Am. Chem. Soc. 106 (1984) 809–810.30 [6] R. Bogseth, E. Edgcomb, C. M. Jones, E. K. Chess, P. Hu, Acetonitrile adduct formation as a sensitive means for simple alcohol detection by31 LC-MS, J. Am. Soc. Mass Spectrom. 25 (2014) 1987–1990.32 [7] S. Pothoczki, L. Temleitner, L. Pusztai, Structure of neat liquids consisting of (prefectly and nearly) tetrahedral molecules, Chem. Rev. 11533 (2015) 13308–13361.34 [8] H. Bertagnolli, P. Chieux, M. D. Zeidler, Neutron-diffraction study of liquid acetonitrile I. CD3C14N., Mol. Phys. 32 (1976) 759–773.35 [9] H. Bertagnolli, P. Chieux, M. D. Zeidler, Neutron-diffraction study of liquid acetonitrile II. CD3C15N., Mol. Phys. 32 (1976) 1731–1736.36 [10] H. Bertagnolli, M. D. Zeidler, Molecular pair-correlation function of liquid acetonitrile from x-ray and neutron-diffraction studies., Mol. Phys.37 35 (1978) 177–192.38 [11] T. Radnai, S. Itoh, H. Ohtaki, Liquid structure of N,N-dimethylformamide, acetonitrile and their 1:1 molar mixture., Bull. Chem. Soc. Jpn.39 61 (1988) 3845–3852.40 [12] S. Pothoczki, L. Pusztai, Intermolecular orientations in liquid acetonitrile: New insights based on diffraction measurements and all-atom41 simulations, J. Mol. Liq. 225 (2017) 160–166.42 [13] M. Kowsari, L. Tohidifar, Systematic evaluation and refinement of existing all-atom force fields for the simulation of liquid acetonitrile, J.43 Comput. Chem. 39 (2018) 1843–1853.44 [14] V. A. Koverga, O. M. Korsun, O. Kalugin, B. Marekha, A. Idrissi, A new potential model for acetonitrile: Insight into the local structure45 organization, J. Mol. Liquid. 233 (2017) 251–261.46 [15] M. Orhan, Dielectric and transport properties of acetonitrile at varying temperatures: a molecular dynamics study, Bull. Korean Chem. Soc.47 35 (2014) 1469–1478.48 [16] M. Albertí, A. Amat, F. De Angelis, F. Pirani, A model potential for acetonitrile: from small clusters to liquid, J. Phys. Chem. B 117 (2013)49 7065–7076.50 [17] A. M. Nikitin, A. Lyubartsev, New six-site acetonitrile model for simulations of liquid acetonitrile and its aqueous mixtures, J. Comput. Chem.51 28 (2007) 2020–2026.52 [18] P. J. Gee, W. van Gusteren, Acetonitrile revisited: a molecular dynamics study of the liquid phase, Mol. Phys. 104 (2006) 477–483.53 :Accepted manuscript Page 16 of 18 Simulation of acetonitrile in liquid, gas and solid phases [19] E. Guàrdia, R. Pinzón, J. Casulleras, M. Orozco, F. Luque, Comparison of different three-site interaction potentials for liquid acetonitrile,1 Mol. Simul. 26 (2001) 287–306.2 [20] X. Grabuleda, C. Jaime, P. A. Kollman, Molecular dynamics simulation studies of liquid acetonitrile: New six-site model, J. Comput. Chem.3 21 (2000) 901–908.4 [21] E. Cabaleiro-Lago, M. Ríos, Intermolecular potential for acetonitrile based on ab initio calculations, Mol. Phys. 96 (1999) 309–321.5 [22] J. Siebers, U. Buck, T. Beu, Calculation of structures and vibrational spectra of acetonitrile clusters, Chem. Phys. 239 (1998) 549–560.6 [23] W. L. Jorgensen, J. M. Briggs, Monte carlo simulations of liquid acetonitrile with a three-site model, Mol. Phys. 63 (1988) 547–558.7 [24] D. Edwards, P. A. Madden, I. McDonald, A computer simulation study of the dielectric properties of a model of methyl cyanide, Mol. Phys.8 51 (1984) 1141–1161.9 [25] H. Böhm, I. McDonald, P. A. Madden, An effective pair potential for liquid acetonitrile, Mol. Phys. 49 (1983) 347–360.10 [26] J. Chen, P. Sit, Ab initio study of the structural properties of acetonitrile-water mixtures, Chem. Phys. 457 (2015) 87–97.11 [27] R. Dias, C. da Costa, T. Manhabosco, A. de Oliveira, M. Matos, J. Soares, R. J. Batista, Ab initio molecular dynamics simulation of methanol12 and acetonitrile:the effect of van der Waals interactions, Chem.Phys. Lett. 714 (2019) 172–177.13 [28] M. Barrow, 𝛼acetonitrile at 215 K, Acta Cryst. B37 (1981) 2239–2242.14 [29] O. K. Antson, K. J. Till, N. H. Andersen, Neutron powder diffraction study of deuterated 𝛽-acetonitrile, Acta Cryst. B43 (1987) 296–301.15 [30] R. Enjalbert, J. Galy, CH3CN: X-ray structural investigation of a unique single crystal. 𝛽→𝛼a phase transition and crystal structure, Acta16 Cryst. B58 (2002) 1005–1010.17 [31] M. Hloucha, A. K. Sum, S. I. Sandler, Computer simulation of acetontrile and methanol with ab initio-based pair potentials, J. Chem. Phys.18 113 (2000) 5401–5406.19 [32] H. Saint-Martin, J. Hernández-Cobos, M. I. Bernal-Uruchurtu, I. Ortega-Blake, H. J. C. Berendsen, A mobile charge densities in harmonic20 oscillators (MCDHO) molecular model for numerical simulations: The water–water interaction, J. Chem. Phys. 113 (2000) 10899–10912.21 [33] H. Saint-Martin, J. Hernández-Cobos, I. Ortega-Blake, Water models based on a single potential energy surface and different molecular22 degrees of freedom, J. Chem. Phys. 122 (2005) 224509(1)–224509(12).23 [34] M. Valdéz-González, H. Saint-Martin, J. Hernández-Cobos, R. Ayala, E. Sánchez Marcos, I. Ortega-Blake, Liquid methanol Monte Carlo24 simulations with a refined potential which includes polarizability, nonadditivity, and intramolecular relaxation, J. Chem. Phys. 127 (2007)25 224507(1)–224507(14).26 [35] E. Galbis, J. Hernández-Cobos, R. Pappalardo, E. Sanchez Marcos, Collecting high-order interactions in an effective pairwise intermolecular27 potential using the hydrated ion concept: The hydration of Cf(III), J. Chem. Phys. 140 (2014) 214104–(1)–(11).28 [36] M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G. A.29 Petersson, H. Nakatsuji, M. Caricato, X. Li, H. P. Hratchian, A. F. Izmaylov, J. Bloino, G. Zheng, J. L. Sonnenberg, M. Hada, M. Ehara,30 K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, J. A. Montgomery, Jr., J. E. Peralta,31 F. Ogliaro, M. Bearpark, J. J. Heyd, E. Brothers, K. N. Kudin, V. N. Staroverov, R. Kobayashi, J. Normand, K. Raghavachari, A. Rendell,32 J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, N. Rega, J. M. Millam, M. Klene, J. E. Knox, J. B. Cross, V. Bakken, C. Adamo, J. Jaramillo,33 R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, R. L. Martin, K. Morokuma, V. G. Zakrzewski,34 G. A. Voth, P. Salvador, J. J. Dannenberg, S. Dapprich, A. D. Daniels, Farkas, J. B. Foresman, J. V. Ortiz, J. Cioslowski, D. J. Fox, Gaussian0935 Revision D.01, 2009. Gaussian Inc. Wallingford CT.36 [37] E. Cabaleiro-Lago, J. Hermida-Ramón, A. Peña Gallego, E. Martínez-Nuñez, A. Fernández-Ramos, Intermolecular interactions and37 cooperative effects in acetonitrile clusters. an based ab initio molecular orbital study, J. Mol. Struct. (Theochem) 498 (2000) 21–28.38 [38] R. A. Mata, B. Costa Cabral, Structural, energetic, and electronic properties of (CH3CN)2–8 clusters by density functional theory, J. Mol.39 Struct. (Theochem) 673 (2004) 155–164.40 [39] K. Remya, C. H. Suresh, Cooperativity and cluster growth patterns in acetonitrile: A dft study., J. Comput. Chem. 35 (2014) 910–922.41 [40] S. N. Ghosh, R. Trambarulo, W. Gordy, Electric dipole moments of several molecules from the stark effect, J. Chem. Phys. 21 (1953) 308–310.42 [41] C. C. Costain, Determination of molecular structures from ground state rotational constants, J. Chem. Phys. 29 (1958) 864–874.43 [42] A. D. Buckingham, R. E. Raab, The dielectric polarization of some imperfect polar gases, J. Chem. Soc. (1961) 5511–5523.44 [43] T. B. Freedman, E. R. Nixon, Matrix isolation studies of methyl cyanide in solid argon and methyl isocyanide, Spectrochim. Acta 28A (1973)45 1375–1391.46 [44] T. A. Renner, M. Blander, A study of dimerization in acetonitrile vapor by measurement of thermal conductivity, J. Phys. Chem. 81 (1977)47 857–861.48 [45] I.T. Todorov,W. Smith,K.Trachenko, M. T. Dove, DL_POLY_3:Newdimensions in molecular dynamics simulations via massiveparallelism,49 J. Mat. Chem. 16 (2006) 1911–1918.50 [46] D. Frenkel, B. Smit, Understanding Molecular Simulation. from Algorithms to Applications, 2nd ed., Academic Press, San Diego,USA, 2002.51 [47] I.-C. Yeh, and Gerhard Hummer, System-Size Dependence of Diffusion Coefficients and Viscosities from Molecular Dynamics Simulations52 with Periodic Boundary Conditions, J. Phys. Chem. B 108 (2004) 15873–15879.53 [48] M. P. Allen, D. J. Tildesley, Computer Simulation of Liquids, Oxford University Press, Oxford, 1987.54 [49] Y. Zhang, A. Otani, E. J. Maginn, Reliable viscosity calculation from equilibrium molecular dynamics simulations: A time decomposition55 method, J. Chem. Theory Comput. 11 (2015) 3537–3546.56 [50] W. E. Putnam, D. M. McEachern, J. E. Kilpatrick, Entropy and related thermodynamic properties of acetonitrile (methyl cyanide), J. Chem.57 Phys. 42 (1965) 749–755.58 [51] S. Pothoczki, L. Temleitner, L. Pusztai, Detailed intermolecular structure of molecular liquids containing slightly distorted tetrahedral59 molecules with c3v symmetry: Chloroform, bromoform and methyl-iodide, J. Chem. Phys. 134 (2011) 044521.60 [52] C. Caleman, P. J. van Maaren, M. Hong, J. S. Hub, L. Costa, D. D. van der Spoel, Force field benchmark of organic liquids: density, enthalpy of61 vaporization, heat capacities, surface tension, isothermal compressibility, volumetric expansion coefficient, and dielectric constant, J. Chem.62 Theory Comput. 8 (2012) 61–74.63 :Accepted manuscript Page 17 of 18 Simulation of acetonitrile in liquid, gas and solid phases [53] M. Kowsari, L. Tohidifar, Tracing dynamics, self-diffusion, and nanoscale structural heterogeneity of pure and binary mixtures of ionic liquid1 1-hexyl-2,3-dimethylimidazolium bis(fluorosulfonyl)imide with acetonitrile: Insights from molecular dynamics simulations, J. Phys. Chem.2 B 120 (2016) 10824–10838.3 [54] X.-W. An, M. Mänsson, Enthalpies of combustion and formation of acetonitrile, J. Chem. Thermodyn. 15 (1983) 287–293.4 [55] DDBST GmbH webpage, Dortmund data bank, 2019. Http://www.ddbst.com/en/EED/PCP/DEN-C3.php.5 [56] R. W. Gallant, Physical properties of hydrocarbons .36. nitriles, Hydrocarbons Processing 48 (1969) 135.6 [57] Y. Marcus, The Properties of Solvents, Wiley, New York, 1998.7 [58] R. Rivelino, B. Costa Cabral, K. Coutinho, S. Canuto, Electronic polarization in liquid acetontrile: A sequential Monte Carlo/quantum8 mechanics investigation, Chem. Phys. Lett. 407 (2005) 13–17.9 [59] Y. Fujita, T. Ohba, S. Ikawa, Measurements of far-infrared optical constants of liquids with a synchrotron radiation source, Can. J. Chem. 6910 (1991) 1745–1751.11 [60] P. Yuan, M. Schwartz, Molecular reorientation in acetonitrile. a comparison of diffusion coefficients from raman bandshapes and nuclear12 magnetic resonance relaxation times, J. Chem. Soc. Farady Trans. 86 (1990) 593–596.13 [61] A. Sugitani, S. Ikawa, S. Konaka, Effect of temperature on the infrared band shapes and reorientational and vibrational relaxation of liquid14 acetonitrile, Chem. Phys. 142 (1990) 423–430.15 [62] H. Shekaari, M. T. Zafarani-Moattar, S. Mirheydari, Density, viscosity, speed of sound, and refractive index of a ternary solution of aspirin,16 1-butyl-3-methylimidazolium bromide, and acetonitrile at different temperatures t= (288.15 to 318.15) k, J. Chem. Eng. Data 60 (2015)17 1572–1583.18 [63] S. Le Roux, V. Petkov, ISAACS – interactive structure analysis of amorphous and crystalline systems, J. Appl. Cryst. 43 (2010) 181–185.19 [64] W. L. Jorgensen, D. S. Maxwell, J. Tirado-Rives, Development and testing of the opls all-atom force field on conformational energetics and20 properties of organic liquids, J.Am. Chem. Soc. 118 (1996) 11225–11236.21 :Accepted manuscript Page 18 of 18