Force field parameterization from the Hirshfeld molecular electronic

Subscriber access provided by Kaohsiung Medical University

Condensed Matter, Interfaces, and Materials

Force field parameterization from the Hirshfeld molecular electronic density Alexander Pérez de la Luz, Jorge Alberto Aguilar-Pineda, Jose Guillermo Mendez-Bemúdez, and José Alejandre J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.8b00554 • Publication Date (Web): 02 Oct 2018 Downloaded from http://pubs.acs.org on October 7, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Force field parameterization from the Hirshfeld molecular electronic density

September 4, 2018 Alexander Pérez de la Luza , Jorge Alberto Aguilar-Pinedaa , José Guillermo Méndez-Berm´ udezb and José Alejandrea∗ a Departamento de Qu´ımica, Universidad Autónoma Metropolitana-Iztapalapa. Av. San Rafael Atlixco 186, Col. Vicentina, 09340, Ciudad de México, México. b Centro Universitario de los Valles, Universidad de Guadalajara, Carretera Guadalajara-Ameca km. 45.4, Ameca 46600, Jalisco, México.

∗

Author to whom correspondence should be addressed. E-mail: [email protected]

ACS Paragon Plus Environment 1

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Abstract The Hirshfeld charges are linearly increased to reproduce the experimental dielectric constant of 10 polar solvents having values between 13 (pyridine) and 182 (N-methylformamide). The OPLS/AA force field is used to obtain the new parameters. The surface tension and liquid density are also target properties to determine the new nonbonding parameters. The charge scaling factor is between 1.2 and 1.3. In addition, properties that were not used in the parameterization procedure, such as the heat of vaporization, self-diffusion coefficient, shear viscosity, isothermal compressibility and volumetric expansion coefficient are obtained. Binary mixtures of amide/water and amide/amide are also studied. The original parameters of OPLS/AA, CGenFF and GAFF force fields are evaluated. The TIP4P/ε force field is used to simulate water. The results from this work with the new parameters, for both pure components and binary mixtures, are in better agreement with experimental data than those obtained with the original values for most of the calculated properties. The maximum density of N-methylformamide in aqueous solutions is correctly predicted only with the new parameters. The high value of the dielectric constant of acetamide, formamide and N-methylformamide is discussed in terms of the chain formation from the hydrogen bond interactions.


Page 2 of 35

Page 3 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1

Introduction

During the last two decades there has been a great interest to develop force fields able to reproduce accurately experimental thermodynamic, transport and structural properties of fluids and solids over a wide range of temperatures and pressures.1, 2, 3, 4 In simulations that include all the atoms, the intramolecular parameters for bond distances, bending angles, torsional barriers and partial atomic charges are generally obtained from quantum mechanical calculations of isolated molecules. The OPLS/AA,1, 3 CGenFF,5 GAFF2 models use Monte Carlo or molecular dynamics simulations to obtaining the van der Waals parameters in the liquid phase to reproduce the density and heat of vaporization. The TraPPE4 and NERD6 models use the coexisting vapor/liquid desities and critical properties. Most of the parameterization procedures are focused in obtaning interaction parameters for pure components of small molecules and in some cases, particularly TraPPE, also binary mixtures. The parameters are transferred to build new molecules and to study other thermodynamic conditions. A benchmark of the OPLS/AA, CGenFF and GAFF force fields for 146 organic liquids as pure components was reported by Caleman et al.7 in 2012 and by Fischer et al.8 in 2015 who used the Ewald sums to calculate the Lennard-Jones interactions. The calculated properties were liquid density, enthalpy of vaporization, heat capacity, surface tension, isothermal compressibility, volumetric expansion coefficient and the dielectric constant. Most of the calculated properties agreed well with the experimental data, but the surface tension and the dielectric constant were systematically different. The values of the dielectric constant were around half of the experimental value. Zubillaga et al.9 and Fischer et al.8 found that long-range interactions from the Lennard-Jones potential affect mainly the surface tension but not the dielectric constant. In 2018, N´ un ˜ez-Rojas et al.10 performed a benchmark of the TraPPE UA force field4 for 41 polar liquids to determine its ability to reproduce the experimental values of the seven properties reported by Caleman et al. The largest difference between the calculated and experimental values was found for the dielectric constant which was underestimated by around 40%. The use of dielectric constant as target property to determine the molecular dipole moment, that involves charges and geometry, has been applied with succes to parameterize water,11, 12 dichloromethane,13 methanol with three14 and four sites,15 formamide (FM),16 acetamide (ACE)17 and formic acid.18 In 2015, Salas et al.14 proposed a systematic procedure to develop force fields for polar liquids. The procedure relates molecular parameters with macroscopic properties. The new non-bonding parameters in this work are obtained using the Salas et al. method14 while the intramolecular parameters are taken from a known force field, such as OPLS/AA, and they are not modified durACS Paragon Plus Environment 3


ing the parameterization procedure. The geometry and charges are used to calculate the molecular dipole moment. The dielectric constant was used to determine the optimum molecular dipole moment by scaling linearly only the partial charges of all the atoms, keeping constant all the LJ parameters. The new charges are fixed and the surface tension at the liquid-vapor interface and the liquid density were used to obtain the Lennard-Jones parameters of energy, εLJ , and size, σLJ , respectively, by scaling linearly the values from the original force field. The Salas method was applied using, as starting point, the parameters of OPLS/AA, CGenFF, GAFF and TraPPE UA force fields for molecules having small dielectric constant.14 We were able to reproduce the three target properties but the selfdiffusion coefficient was systematically lower for most of the studied systems compared with the results obtained with the original parameters and with those of the experimental data. The charge scaling method was applied to liquid FM16 using the OPLS/AA parameters but it was not possible to increase the calculated dielectric constant from 50, that was predicted with the original parameters, to the experimental value of 110. Pérez de la Luz et al.16 showed that by using the solvation model density,19 SMD, in the electronic structure calculations with the Mulliken population analysis20 to obtain the atomic charges it was possible to obtain the experimental value after applying the Salas et al. method to determine the LJ parameters. The use of NBO,21, 22, 23, 24 CHELPG25 and Merz-Kollman,26, 27 MK, schemes with implicit solvent using SMD to obtain the partial charges on FM gave systems in the solid state in simulations at ambient conditions. The charges from MK scheme for the more electronegative atoms, oxygen and nitrogen, were around 1.5 and 2.5 times larger than those obtained from the Mulliken analysis. The new parameters were able to reproduce the experimental values of the target properties as a function of temperature and it showed the importance of the dielectric constant in predicting the correct miscibility of two polar liquids on binary systems. That work showed that the charge distribution has an important effect on the self-diffusion coefficient and other physicochemical properties. Aguilar-Pineda et al.17 parametrized acetamide using the Salas et al. method and they were able to reproduce the increase of dielectric constant on acetamide-water mixtures only with the new OPLS/AA parameters. We applied in this work the scaling method of Salas et al. to N-methylformamide (NMF), which has a dielectric constant of 182, at room conditions using the Mulliken population analysis to obtain the charges. The self-diffusion coefficient was close to zero and it was not possible to calculate the dielectric constant because the total dipole moment of the system did not converge. One of the main problems in developing force fields that include all the atoms is the assignment of the partial charges in every atom.28 Apart from the methods ACS Paragon Plus Environment 4

Page 4 of 35

Page 5 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


mentioned above to obtain their values, it is also possible to use the partition of molecular electron density of atoms in molecules, AIM, such as the Bader29 and Hirshfeld30 methods. The atomic charges are not an observable of the system, ie, there is not an operator associated with that property and therefore, it is not possible to obtain a unique value using those approximations. In the Hirshfeld scheme, the molecular electron density, ρmol (r), at point r obtained by any computation of electronic structure, is partitioned into charge density contributions of the different atoms bound in the molecule ρb.a. i (r),

ρb.a i (r) =

ρat i (r) mol ρ (r) pro ρ (r)

(1)

here, ρat i (r) are the electronically optimized, spherically averaged densities of isolated pro atoms obtained at the same level of theory as the molecule P at under study and ρ (r) is the electronic density of the promolecule defined as i ρi (r), in this way the charge for the atom i, qi , is calculated as, Z qi = Zi − ρb.a. (2) i (r)dr where Zi is the nuclear charge of the atom i. The partial charges from quantum calculations are generally scaled to be used in condensed matter simulations as it is done in the CMx (x=1,2,3,4,5) charge models.31, 32, 33, 34, 35, 36 In particular, the charges from the Hirshfeld scheme are modified by the CM5 procedure and additionaly they are scaled empirically by a factor of 1.20 to obtain reliable results of density and free energy of hydation37 of several polar liquids. Cole et al.38 proposed a method to obtain environment-specific partial atomic charges and LJ parameters using a linear-scaling density functional theory and atomsin-molecule electron density partitioning different than the Hirshfeld scheme. Virtual sites were added to reproduce the dipole and quadrupole moments of the AIM atomic electronic densities in polar moleculas such as amines. The nonbonding parameters were used with the OPLS/AA force field to obtain the density, heat of vaporization and free energy of hydration of several polar liquids. Hereafter, the force field used by Cole et al. is denoted as OPLS16. In that work, only the LJ energy parameters of some atoms were made transferable. Macchiagodena et. al39 used the CM5 charges for amides with additional virtual sites to model the lone pairs on carbonyl oxigen atom to favor the hydrogen bond interactions. The low value of self-diffusion coefficient obtained with the Salas et al. method14 in our previous calculations might be related with the large value of the partial charges for the more electronegative and electropositive atoms obtained after the scaling ACS Paragon Plus Environment 5


procedure where the attraction between atoms with opposite charge is enhanced. Therefore, it is justified to look for other schemes to calculate the atomic charges. The main goal of this work is to determine if the molecular electronic density distribution from the Hirshfeld scheme allows obtaning the partial atomic charges of polar liquids to be used in condensed phase simulations. In addition, to analyze if it is possible to apply the Salas et al. method to obtain the intermolecular force field parameters that reproduce not only the dielectric constant, surface tension and liquid density but also other thermodynamic and transport properties for 10 polar liquids. It is well known that the Hirshfeld scheme gives very small values for atomic charges compared with Mulliken, NBO, CHELPG and MK methods. The Hirshfeld charges were scaled by a factor between 1.2 and 1.3 to reproduce the experimental dielectric constant. This parameterization procedure is different from the work of Doda et al.37 and Cole et al.38 Virtual sites are not added to any atom in the molecules. The 10 polar liquids studied are shown in Fig. 1. Apart from the intrinsic interest of parameterization of those liquids, some of them are also important in several applications. Propylene carbonate is widely used as a solvent in electrolyte solutions of lithium ion batteries. The molecules that contain the amide group are important to study hydrogen bonding interactions in order to understand the folding process of proteins. The work is organized as follows: The force field and computational details are given in Section 2. The results for pure components and binary mixtures are discussed in Section 3. Finally, the Conclusions and References are given.

2

Force Field and Computational Details

The OPLS/AA force field1 is used to obtain the new parameters for the molecules shown in Fig. 1. The functional form for the nonbondig interactions contains Lennard-Jones and Coulombic terms, ( " ) 6 # 12 σij σij qi qj V (rij ) = 4εij − + fij (3) rij rij 4πε0 rij where rij is the distance between atoms i and j, qi is the partial charge on atom i, ε0 is the vacuum permittivity, σij and εij are the effective atom diameter and interaction strength, respectively. The factor fij = 1.0, except for intramolecular interactions between atom pairs separated by three bonds1, 3, 40 where the value is 0.5. The cross-interactions are calculated using the Lorentz-Berthelot mixing rules: σij = (σii + σjj )/2 and εij = (εii εjj )1/2 . The original intramolecular and nonbonding parameters for the OPLS/AA, GAFF and CGenFF are taken from the web page http : //virtualchemistry.org.41 The paACS Paragon Plus Environment 6

Page 6 of 35

Page 7 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 1: Polar systems studied in this work. The following color code is used: White for hydrogens, gray for carbons, blue for nitrogens and red for oxygens. rameters for the OPLS16 model are taken from the Supporting Information of the work of Cole et al.38 The optimum values of charges and LJ parameters in this work are obtained using the Salas et al.14 procedure where the intramolecular parameters of the OPLS/AA force field are left unchanged. The quantum mechanical computations were performed using the Gaussian09 set of programs42 with the M062X functional43 and the 6-311++g** basis set.44, 45 The SMD was used as it is implemented in Gaussian09. The partial atomic charges were calculated using the standard Hirshfeld30 method. The Hirshfeld charges were scaled linearly to reproduce the experimental value of the dielectric constant of solvents using molecular dynamics simulations. The MD calculations for liquids are performed at temperatures where experimental properties are available at 1 bar. The results obtained with the new parameters are denoted as NewFFP. A table with the original OPLS/AA and new non-bonding parameters for every atom in a molecule is included as a Supporting Information.



The molecular dynamics simulations were carried out using the GROMACS 4.5.4 package46 with the leap-frog algorithm and a time step of 2 fs to solve the equations of motion using periodic boundary conditions in the three directions. In all the simulations, the bond distances are kept rigid by using the LINCS alogrithm.47 The electrostatic interactions are obtained with the Particle Mesh Ewald approach48 with a tolerance of 10−6 for the real space contribution and with a grid spacing of 1.2 ˚ A and spline interpolation of order 4 for the reciprocal space. The NPT ensemble (constant number of molecules, temperature and pressure) with isotropic fluctuations of volume is used to calculate the properties in the liquid phase. The temperature was coupled to the Nosé-Hoover thermostat with a parameter τT =0.5 ps while the pressure was coupled to the Parrinello-Rahman barostat with a coupling parameter τP =1.0 ps. Those values are adequate to sample properly the volume fluctuations.10 These simulations involved typically 500 molecules. The real part of the Ewald summation and the LJ interactions were truncated at 12 ˚ A. Long range corrections for the LJ energy and pressure were added. The dielectric constant was obtained using the dipole moment fluctuations of the system.49, 50 The liquid properties were calculated for at least 150 ns after an equilibrition period of 20 ns. The NVT ensemble (constant number of molecules, volume and temperature) is used to obtain the surface tension at the liquid-vapor interface using a slab geometry. Initially a liquid containing 2500 molecules is equibrated at the temperature of interest and 1 bar using the NPT ensemble and then the simulation cell is elongated in the z direction to decrease the density of the system in order to put it in the liquid-vapor region. At equilibrium, the liquid coexists with a vapor in the longest direction of the simulation box. Typical values of the simulation cell dimensions were Lx = Ly =52 ˚ A with Lz = 3Lx , where z being the normal direction to the liquid-vapor interface. The cutoff distance was set to 25 ˚ A to avoid fine size effects on truncation distance51, 9 52, 53 and on the cross sectional area. Long range corrections were not added for the LJ energy and pressure. The surface tension was calculated using the mechanical definition with the diagonal components of the pressure tensor.51, 9 The equilibration period for the interfacial simulations was 10 ns and the results for the average properties were obtained over an additional 10 ns trajectory.

3

Results

Apart from the target properties, the heat of vaporization, isothermal compressibility, volumetric expansion coefficient, self-diffusion coefficient and shear-viscosity are calculated for pure components. We used the relative error, |Xexp − XM D | ∗ 100/Xexp , to evaluate the differences between calculated properties with the new parameters from ACS Paragon Plus Environment 8

Page 8 of 35

Page 9 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


OPLS/AA with those from experimental data. The same analysis was made for calculations with the original parameters of different force fields (OPLS/AA,1 GAFF,2 CGenFF5 and OPLS1638 ). The chain formation of ACE, FM and NMF is discussed in terms of the hydrogen bond interactions. The liquid density and dielectric constant are also calculated for selective amide/water and amide/amide binary mixtures. The equations used to calculate the properties are given as Supporting Information.

3.1 3.1.1

Simulations of Pure Liquids Dielectric Constant, Surface Tension and Liquid Density

The dielectric constant is obtained using the dipole moment fluctuations of the system.50 The correlation between the calculated results using different force fields with those from the experimental data are shown in Figure S1-A as Supporting Information and the numerical values are given in Table S1. The results from this work are in excellent agreement with experimental data but those from other force fields with the original parameters are systematically lower for most of the systems. For the OPLS16 model, the largest deviation is for acetamide where the value is around half of the experimental data. The average relative error with the new parameters is 3.4 while the largest value is 52.6 for the OPLS/AA model, see Table 1. The largest differences in dielectric constant for all the force fields with the original parameters are found for ACE, FM and NMF. The values reported by Caleman et al.7 for NMF at 298.15 K and 1 bar using GAFF and OPLS/AA were 13.3 ± 0.3 and 19.2 ± 0.7, respectively, compared with the experimental value of 182. Those results are shown in Figure S1-A with an open square and a open diamond. From the coordinates provided in the web page http : //virtualchemistry.org is found that the oxygen atom is in a cis conformation (zero degrees for dihedral angle OCNH) with the hydrogen bonded to a nitrogen atom. From NMR experiments63 is found that NMF is a mixture of 95% of molecules in the trans conformation (180 degrees of dihedral angle OCNH) for those atoms. The molecular dipole moment with the charges of GAFF model increases from 4.23 D to 4.69 D when the conformation is changed from cis to trans. The results obtained in this work with the original parameters for OPLS/AA and GAFF force fields but with the trans conformation were 68 and 191.9, respectively. Those results highlight the importance of the molecular dipole moment, that includes the charge distribution and geometry, to calculate the dielectric constant. As we will see below, the high value of the dielectric constant is related with the formation of chains due to the hydrogen bond interactions between the NH · · · OC atoms from two molecules, rather than, to the increase of the dipole moment of a single molecule. The results for FM and NMF obtained by Macchiagodena et. al39 using the OPLS/AA model with virtual sites on the oxygen atom of the carbonyl group are shown in Fig. S1-A with open circles and they are in good agreement with experiment. The results for ACS Paragon Plus Environment 9


both molecules using a polarizable model reported by Xie et. al56 are overestimated by 26% for FM and 11% for NMF and they are shown with stars in Fig. S1-A. The surface tension is obtained using the diagonal components of the pressure tensor.51, 9 The correlation between the calculated results using different force fields with those from the experimental data are shown Fig. S1-B as Supporting Information and the numerical values are given in Table S2. The results from this work with the new parameters are in excellent agreement with experimental data and the average relative error for all the molecules is 2.7, see Table 1. The calculated values obtained with the other force fields are scattered and systematically overestimated where the largest average relative error is 26.1 for GAFF. It is surprinsing that the results from all the force fields for acetone agree well with experiment. The largest error for individual molecules is 74% for nitromethane with GAFF. The correlation between the calculated liquid density using different force fields with those from the experimental data are shown Fig. S1-C as Supporting Information and the numerical values are given in Table S3 . The average relative error is 0.3 for results with the new parameters, see Table 1. The liquid density is also a target property to parameterize the other force fields used in this work, therefore, it is expected to have good agreement with experimental data. The largest deviation is for GAFF where the relative error is 4.1. 3.1.2

Heat of Vaporization, Self-Diffusion Coefficient and Shear-Viscosity

The heat of vaporization is obtained using the difference between the entalphy of the vapor and liquid assuming the vapor behaves as an ideal gas. The same procedure has been used to calculate that property to parameterize the original force fields used in this work. The polarizable corrections are not included. In our case the heat of vaporization is not a target property while in the other force fields it is. It is a good property to evaluate the quality of the new parameters to calculate a property that was not included in the parameterization procedure. The correlation between the calculated heat of vaporization with those from the experimental data are shown Fig. S2-A as Supporting Information and the numerical values are given in Table S4 . The largest difference, 22 %, for an individual molecule is for acetamide where the simulation was performed at 364.15 K and 1 bar because it is solid below that temperature. The vapor pressure, the normal pressure to the interface, at 364.15 K was determined in the same simulation performed to calculate the surface tension, then, a simulation in the NPT ensemble was performed in the vapor phase. The result for the entalphy was essentially the same as that obtained with the ideal gas approximation. The average relative error for the heat of vaporization is 4.5 for results with the new parameters, see Table 1, that value is slightly larger than 3.4 calculated for the OPLS/AA. The largest deviation is for GAFF where the average error is 16.3. ACS Paragon Plus Environment 10

Page 10 of 35

Page 11 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The self-diffusion coefficient is obtained using the mean square displacement correlation function through the Einstein equation. The correlation between the calculated results using different force fields with those from the experimental data are shown in Fig. 2S-B as Supporting Information and the numerical values are given in Table S5. The average relative error is 19.2 for results with the new parameters, see Table 1. The agreement with experiment is excellent using the new parameters considering that this property is not used in the parameterization procedure. The largest difference, 28%, is observed for acetone. The largest relative error is 60.6 for the CGenFF. The simulation results are underestimated for all of the force fields. The shear-viscosity is obtained using the correlation function of the non-diagonal components of the pressure tensor. The correlation between the calculated results using different force fields with those from the experimental data are shown Fig. 2SC as Supporting Information and the numerical values are given in Table S6 . The average relative error is 21.8 for results with the new parameters, see Table 1. The formamide and propylene carbonate are the systems where the shear-viscosity is 32% and 43% greater than the experimental values. The best and the worst agreement is for OPLS/AA, 17%, and 52% for GAFF, respectively. 3.1.3

Isothermal Compressibility and Volumetric Expansion Coefficient

The isothermal compressibility is obtained using the volumen fluctuations in NPT simulations. The correlation between the calculated results using different force fields with those from the experimental data are shown Fig. S3-A as Supporting Information and the numerical values are given in Table S7. The average relative error is 11.6 for results with the new parameters, see Table 1. The smallest error, 8.7, is for the OPLS/AA force field and the largest value is 27.5 for GAFF. The results are scattered for the four systems studied with the OPLS16 model where the average relative error is 25.5. The volumetric expansion coefficient is obtained using the termodynamic definition, αT = −(∂ρ/∂T )P /ρ. The density was obtained for 5 temperatures separated by 10 K in NPT simulations at 1 bar. Two values above and two values below of the target temperature, except for acetamide where the temperature was varied in 5 degrees because it is a solid close to the target temperature. The liquid density was fitted to a linear equation and the value of αT was interpolated from those results. The correlation between the calculated values using different force fields with those from the experimental data are shown in Fig. S3-B as Supporting Information and the numerical values are given in Table S8 . The average relative error is 15.5 for results with the new parameters, see Table 1.



Table 1: Relative error, |Xexp − XM D | ∗ 100/Xexp , between calculated and experimental values for the 10 polar liquids studied in this work. The values with green and red colors are for results closest and farthest from experimental values, repectively. Property NEWFF OPLS/AA GAFF CGenFF OPLS16 Dielectric Constant 3.4 52.6 37.6 42.8 28.8 Surface Tension 2.7 10.5 26.1 22.1 20.7 Density 0.29 1.8 4.1 3.2 2.6 Heat of Vaporization 4.5 3.4 16.3 9.9 5.7 Self-Diffusion 19.2 37.4 49.0 60.6 32.9 Shear Viscosity 21.8 16.6 52.2 35.2 14.0 Isothermal Compressibility 11.6 8.7 27.5 19.2 25.5 Volumetric Expansion Coefficient 15.5 6.2 8.0 4.0 13.4 3.1.4

Molecular Self-Assembly Due to Hydrogen Bond Interactions Between Amides

To analyze the formation of chains between amide molecules we calculated the pair distribution function between the oxygen atom with a hydrogen bonded to a nitrogen atom (NH. . . OH). The oxygen and hydrogen atoms within the same molecule are in trans conformation. Those atoms have hydrogen bond interactions with other molecules. The results for ACE, FM and NMF are shown in Fig. 2. The first peak is for intermolecular interactions while the second and third are for intramolecular. In NMF there is only one O-H pair within the same molecule while ACE and FM have two. The strong interaction is for NMF followed by ACE and FM. The reason why acetamide is stronger than formamide is because those molecules favor the formation of rings, that have two hydrogen bonds at around the same distance, rather than chains where only one hydrogen bond is observed. 8

Acetamide Formamide N-Methylformamide

6 g(r)(NH...OH)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 35

4

2

0 0.1

0.2

0.3 r[nm]

0.4

0.5

Figure 2: Pair distribution function for oxygen with hydrogens bonded to a nitrogen atom in ACE, FM and NMF. The first peak is for intermolecular interactions while the second and third peaks are for intramolecular.


Page 13 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The dipole moment of a single molecule is 5.41 D for ACE, 4.89 D for FM and 5.16 D for NMF. These results can not explain why the dielectric constant has values of 67, 109 and 182 for those molecules. Clearly, there is a cooperative effect due to the hydrogen bond interactions. To analyze the formation of chains in those molecules we applied a geometric criteria. Two molecules are part of a chain if the oxygen and hydrogen atoms from two molecules are at a distance slightly greater than that of the first peak position of the pair distribution function. The method of Sevick80 was used to determine the length of the chains. This approach does not give the total number of hydrogen bonds in the whole system but it is useful to have a rough idea about the effect of chain length has on the dielectric constant. The amide systems can be seen as a mixture of clusters or new molecules with dipole moments larger than that of a single molecule. For instance, two clusters of formamide taken from a liquid simulation contain 12 and 13 molecules with dipole moments of 19 D and 35 D, respectively. The first one has a compact structure while the second one is linear.

Figure 3: Typical chains of amides from simulations at 1 bar: A) Acetamide at 364.15 K. B) Formamide at 298.15 K and C). N-methylformamide at 298.15 K. The following colors are used to distinguish every atom: white for hydrogen, red for oxygen, blue for nitrogen and gray for carbon. The hydrogen bond is defined in terms of the oxygen atom of the carbonyl group with a hydrogen bonded to a nitrogen atom. The electrostatic potential is a guide to the eye to follow the hydrogen bonding network. The molecular electrostatic potential was obtained using the Gaussian09 program for a chain extracted from a liquid. The images in Fig. 3 for the three amides are used to show visually the interaction between the NH· · · OC of two different molecules to form hydrogen bonds. The structure of ACE, see Fig. 3-A, and FM, see Fig. 3-B, has an oxygen atom that share two hydrogens bonded to a nitrogen making the cluster more compact. In the case of NMF, see Fig. 3-C, one oxygen from one molecule form one hydrogen bond with other molecule, this makes larger the dipole moment of the chain compared with those of the other two amides. ACS Paragon Plus Environment 13


Figure 4 shows a typical cluster size distribution at 1 bar for liquid ACE at 364 K, FM and NMF at 298.15 K from NPT simulations using 500 molecules for one equilibrated configuration. We observe a large number of molecules that are not connected with hydrogen bonds. The number is larger for NMF and smaller for FM, in agreement with the results of the radial distribution function. The NMF and ACE form clusters containing up to 22 and 11 molecules, respectively. The number of cluster with small size is larger in ACE as it is seen in the inset of Fig. 4. 10 8 6 4 2 0 -2

200

Number of clusters

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 35

150

0

5

10

15

20

25

100 Acetamide Formamide N-methylformamide 50

0

0

5

10 15 Zize of clusters

20

25

˚), FM (1.89 ˚ Figure 4: Cluster size distributions for ACE(1.84 A A) and NMF(1.74 ˚ A). The numbers between parenthesis are the used distances to determine if two molecules are connected by hydrogen bonds. The values are slightly larger that the position of the first maximum in the radial distribuction functions shown in Fig. 2.

3.2

Binary Mixtures Simulations

The new parameters were obtained for pure components. An important test is to use them to study binary mixtures to evaluate if the they are transferable to reproduce experimental properties at different compositions. The TIP4P/ε model is used because reproduces the experimental dielectric constant and several density anomalies.


Page 15 of 35

3.2.1

Liquid Density of Amide/water Solutions

The results for liquid density of mixtures between water/ACE, water/FM and water/NMF are calculated at 298.15 K and 1 bar and they are shown in Figure 5 as a function of molar fraction of the amide, X(Amide). At that temperature, ACE is solid as a pure component but it disolves in water at room conditions. In 2015, Aguilar et al.17 parameterized ACE using the Salas et al. method14 with charges from the Mulliken analysis and LJ parameters from OPLS/AA. The methyl was treated as united atom. The parameters were used to study acetamide in aqueous solutions with the TIP4P/ε model for water. They found that ACE in solution formed rings and chains due to hydrogen bond interactions between the oxygen atoms and the hydrogens bonded to a nitrogen in different molecules. The experimental density of ACE (1160 kg/m3 ) and FM (1128 g/m3 ) at room conditions is greater than that of water (996 kg/m3 ) and therefore that property increases as the amide concentration increases, see Fig. 5-A and Fig. 5-B. All the force fields predict the correct trend but the best agreement is for the new and original parameters using OPLS/AA. The results for NMF are shown in Fig. 5-C. The experimental density for that liquid is 999 kg/m3 . The maximum is reproduced with the new set of parameters. The failure of the other force fields to reproduce the correct behavior with composition might be due to the large deviations of the liquid density for the pure component. 1300

1060

1075 B)

C)

A) 1050

Exp NewFF OPLS/AA GAFF CGenFF

1250

1040

1050

1200

3

ρ [kg/m ]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1030 1150

1025

1020 1100 1010 1000 1050 1000

990

1000 0

10 20 X(ACE)

30

0

20

40 60 X(FM)

80

100

975

0

20

40 60 80 X(NMFM)

100

Figure 5: Liquid density of binary mixtures as a function of amide concentration at 298.15 K and 1 bar. A) Acetamide/water. B) Formamide/water. C) Nmethylformamide/water. All the experimental data are taken from Jelinska et al.81 and Papama et al.82 The simulation error for the three properties is less than the symbol size.



3.2.2

Dielectric constant of amide/water and amide/amide

The results for dielectric constant of ACE/water, FM/water and NMF/water and FM/NMF mixtures are shown in Figure 6. The new parameters predict correctly the increase of that property for the four mixtures. All the force fields with the original parameters predict an incorrect trend for all the systems except the GAFF model applied to ACE/water and FM/NMF mixtures. The results using GAFF are expected because that model predict the dielectric constant in good agreement with experimental data. 100

125 A)

B)

90

100

80 ε

75 70 50

60 50

0

10

250 C)

200 150

20 X(ACE)

30

40

Exp Exp NewFF OPLS/AA GAFF CGenFF

25

0

200

20

40 60 X(FM)

80

100

40 60 X(NMF)

80

100

D)

150

ε

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 35

100

100

50

50 0

0

20

40 60 X(NMF)

80

100

0

0

20

Figure 6: Dielectric constant of binary mixtures as a function of amide concentration at 1 bar. A) Acetamide/water at 298.15 K. B) Formamide/water at 298.15 K. C) N-methylformamide/water at 293.15 K. D) Formamide/N-methylformamide at 303.15 K.83, 84, 85, 58 The superscripts are for references of the experimental data. The simulation error for the three properties is less than the symbol size.

3.2.3

Excess entalphy and volume of mixing

The deviation from an ideal solution can be studied using some excess mixing properties.86 The experimental information of binary mixtures is much more limited compared with that of the pure components. However, for the mixtures studied in this work, results for the excess entalphy of FM/NMF and excess volume of water/NMF are available. The excess enthalpy of mixing was calculated using the enthalpy of the mixture and those of the pure components as, ∆HM ix = HM ix − [X(F M )HF M + X(N M F )HN M F ], where the subscripts are for formamide and N-methylformamide, respectively. The excess volume of water/NMF mixtures was calculated ACS Paragon Plus Environment 16

Page 17 of 35

by changing in the previous equation H by V and FM by water. The calculated values for different force fields and experimental data are shown in Fig. 6-A and Fig. 6-B for ∆HM ix and ∆VM ix , respectively. The best agreement in both properties is for the NewFFP. The GAFF and CGenFF predict the correct negative sign but their values are overestimated while the OPLS/AA model predict a positive change. 3

B)

A) 0

-1

3

∆VMix.[cm /mol]

2

∆HMix.[kJ/mol]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1

-2

0

Exp. NEWFFP OPLS/AA GAFF CGenFF

-3

-1

0

0.2

0.4 0.6 X(NMF)

0.8

1

-4

0

0.2

0.4 0.6 X(NMF)

0.8

1

Figure 7: Results for binary mixtures for different force fields. A) Excess enthalpy of mixing for FM/NMF at 308 K. The experimental data are taken from A. M. Zaichikov.87 B) Excess volume of mixing for water/N at 318 K. The experimental data are taken from J. Zielkiewicz.88 The simulation error is less than the symbol size. Supporting Information This material is available free of charge via internet at http://pubs.acs.org/. Equations of calculated properties, numerical data and figures for dielectric constant, surface tension, liquid density, heat of vaporization, self-diffusion, shear-viscosity, isothermal compressibility and volumetric expansion coefficient. The intramolecular and intermolecular parameters for the 10 polar liquids and the binary mixtures are provided in the format of Gromacs input. The original OPLS/AA and new (NewFFP) non-bonding parameters for every atom in all the molecules are given in Tables. ACKNOWLEDGMENTS We thank to Laboratorio de Supercómputo y Visualización en Paralelo for the time provided to develop this work, to Professor José Luis Gázquez for helpful discussions about the Hirshfeld method and to Conacyt for financial support under the project ACS Paragon Plus Environment 17


257422. APL thanks Conacyt for a PhD scholarship.

4

Conclusions

We investigate a new strategy to obtain the partial charges of atoms to be used in molecular simulations of polar fluids in condensed matter. The total molecular electronic density of atoms from the Hirshfeld partitioning scheme is used. The quantum mechanical calculations are performed using the SMD to include implicitily the solvent effects. The charges obtained in that way give calculated dielectric constants systematically smaller than the experimental values. The initial values of the LJ parameters to perfom the molecular dynamics simulations are taken from the OPLS/AA force field and they are kept constant during the dielectric constant calculations. The Hirshfeld charges are scaled linearly to reproduce the experimental dielectric constant. The increase of the charge values to be used in simulations of the liquid phase can be seen as an effective polarization to take into account the electric field around a molecule. The scaling factors for charges and Lennard-Jones parameters were different for every molecule. The Hirshfeld charges were increased by a factor between 1.2 to 1.3, the εLJ decreased in around 25% and the change in σLJ was around 1.5%. We did not attemp to use the same factor for a homologous group. The scaling factor is similar to the value of 1.2 used for CM5 charges in order to obtain the liquid density and free energy of hydration in polar liquids. Virtual sites were not located in any of the atoms with lone pairs. The scaled Hirshfeld charges are used in conjuntion with the Salas et al.14 method to obtain the LJ parameters through the calculation of the surface tension and liquid density. The new nonbonding parameters improve the results of the original OPLS/AA force field for the three target properties and self diffusion coefficient, see Table 1 for the reative errors. The relative error for the heat of vaporization with the new parameters, that is not a target property in this work, is 4.5 against 3.5 obtained with the original parameters. The viscosity and isothermal compressibility values are slightly better for calculations with the OPLS/AA with the original parameters. The GAFF model gives the largest relative errors for 5 of the 8 calculated properties for the pure components, see Table 1. The dielectric constant, surface tension, selfdiffusion coefficient and shear-viscosty are the properties with larger relative error using the original parameters for the OPLS/AA, GAFF, CGenFF. The large value of the dielectric constant of amide molecules can not be explained in terms of the molecular dipole moment. The strong hydrogen bond interactions between the NH · · · OC atoms from two molecules form chains with different size, see Fig. 2, where the dipole moment of the chain is larger than that of the single molecule and therefore there is a collective effect. The largest chains are for NMF. In all the ACS Paragon Plus Environment 18

Page 18 of 35

Page 19 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


cases, the formation of chains are favored rather than the cyclic structure. The results from this work show that it is possible to obtain the dielectric constant of highly polar fluids by using a scaled charge distribution from the Hirshfeld partitioning scheme. The new set of parameters give results of liquid density and dielectric constant for amide/water and amide/amide mixtures in excellent agreement with experimental data. The maximum density of NMF in aqueous solutions is also reproduced using those parameters. We did not attempt to transfer parameters from one group of atoms in a molecule to build new structures. The transferability of non-bonding parameters to other molecules and thermodynamics conditions (temperature and pressure) requires additional efforts. The new Lennard-Jones parameters, that come from the OPLS/AA force field, might have the same lack of success than those from the original force field.



References [1] Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids. J. Am. Chem. Soc. 1996, 118, 11225-11236. [2] Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. Development and testing of a general amber force field. J. Comput. Chem. 2004, 25, 11571174. [3] Jorgensen, W. L.; Swenson, C. J. Optimized Intermolecular Potential Functions for Amides and Peptides. Structure and Properties of Liquid Amides. J. Am. Chem. Soc. 1985, 107, 569-578. [4] TraPPE: Transferable Potentials for Phase Equilibria Force Field. http://chemsiepmann.oit.umn.edu/siepmann/trappe/index.html (accessed July 15, 2018). [5] Vanommeslaeghe, K.; Hatcher, E.; Acharya, C.; Kundu, S.; Zhong, S.; Shim, J.; Darian, E.; Guvench, O.; Lopes, P.; Vorobyov, I.; Mackerell Jr, A. D. CHARMM General Force Field: A Force field for Drug-Like Molecules Compatible with the CHARMM All-Atom Additive Biological Force Field. J. Comput. Chem. 2010, 31, 671-690. [6] Nath, S. K.; Escobedo, F. A.; de Pablo, J. J. On the Simulation of vapor-liquid equilibria for alkanes. J. Chem. Phys. 1998, 108, 9905. [7] Caleman, C.; van Maaren, P. J.; Hong, M.; Hub, J. S.; Costa, L. T.; van der Spoel, D. Force Field Benchmark of Organic Liquids: Density, Enthalpy of Vaporization, Heat Capacities, Surface Tension, Isothermal Compressibility, Volumetric Expansion Coefficient, and Dielectric Constant. J. Chem. Theory Comput. 2012, 8, 61-74. [8] Fischer, N. M.; van Maaren, P. J.; Ditz, J. C.; Yildirim, A.; van der Spoel D. Properties of Organic Liquids when Simulated with Long-Range Lennard-Jones Interactions. J. Chem. Theory Comput. 2015, 11, 2938-2944. [9] Zubillaga, R. A.; Labastida, A.; Cruz, B.; Mart´ınez, J. C.; Sánchez, E.; Alejandre, J. Surface Tension of Organic Liquids Using the OPLS/AA Force Field. J. Chem. Theory Comput. 2013, 9, 1611-1615. [10] N´ un ˜ez-Rojas, E.; Aguilar-Pineda, J. A.; Pérez de la Luz, A.; de Jes´ us González, E. N.; Alejandre J. Force Field Benchmark of the TraPPE UA for Polar Liquids: Density, Heat of Vaporization, Dielectric Constant, Surface Tension, Volumetric


Page 20 of 35

Page 21 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Expansion Coefficient, and Isothermal Compressibility. J. Phys. Chem. B. 2018, 122, 1669-1678. [11] Alejandre, J.; Chapela, G. A.; Saint-Mart´ın, H.; Mendoza, N. A non-polarizable model of water that yields the dielectric constant and the density anomalies of the liquid: TIP4Q. Phys. Chem. Chem. Phys. 2011, 44, 19728-19740. [12] Fuentes-Azcatl, R.; Alejandre, J. Non-Polarizable Force Field of Water Based on the Dielectric Constant: TIP4P/ε. J. Phys. Chem. B. 2014, 118, 1263-1272. [13] Fennell, C. J.; Li, L.; Dill, K. A. Simple Liquid Models with Corrected Dielectric Constants. J. Phys. Chem. B. 2012, 116, 6936-6944. [14] Salas, F. J.; Méndez-Maldonado, G. A.; N´ un ˜ez-Rojas, E.; Aguilar-Pineda, G. E.; Dom´ınguez, H.; Alejandre, J. Systematic Procedure To Parametrize Force Fields for Molecular Fluids. J. Chem. Theory Comput. 2015, 11, 683-693. [15] Mart´ınez-Jiménez, M.; Saint-Mart´ın, H. A Four-Site Molecular Model for Simulations of Liquid Methanol and Water-Methanol Mixtures: MeOH-4P. J. Chem. Theory Comput. 2018, 14, 2526-2537. [16] Pérez de la Luz, A.; Méndez-Maldonado, G. A.; N´ un ˜ez-Rojas, E.; Bresme, F.; Alejandre, J. A New Force Field of Formamide and the Effect of the Dielectric Constant on Miscibility. J. Chem. Theory Comput. 2015, 11, 2792-2800. [17] Aguilar-Pineda J. A.; Méndez-Maldonado, G. A.; N´ un ˜ez-Rojas, E.; Alejandre J. Parametrisation of a force field of acetamide for simulations of the liquid phase. Mol. Phys. 2015, 113, 2716-2724. [18] Salas, F. J.; N´ un ˜ez Rojas, E.; Alejandre, J. Stability of formic acid/pyridine and isonicotinamide/formamide cocrystals by molecular dynamics simulations. Theor Chem Acc. 2017, 136, 17. [19] Marenich, A. V.; Cramer, C. J.; Truhlar, D. G. Universal solvation model based on solute electron density and on a continuum model of the solvent defined by the bulk dielectric constant and atomic surface tensions. J. Phys. Chem. B. 2009, 113, 6378-6396. [20] Mulliken, R. S. Electronic Population Analysis on LCAO-MO Molecular Wave Functions. I. J. Chem. Phys. 1955, 23, 1833-1840. [21] Foster, J. P.; Weinhold F. Natural hybrid orbitals. J. Am. Chem. Soc. 1980, 102, 7211-7218.



[22] Reed, A. E.; Weinhold, F. Natural bond orbital analysis of near-Hartree-Fock water dimer. J. Chem. Phys. 1983, 78, 4066-4073. [23] Reed, A. E.; Weinstock, R. B.; Weinhold, F. Natural-population analysis. J. Chem. Phys. 1985, 83, 735-746. [24] Reed, A. E.; Weinhold, F. Natural Localized Molecular Orbitals. J. Chem. Phys. 1985, 83, 1736-1740. [25] Breneman, C. M.; Wiberg, K. B. Determining atom-centered monopoles from molecular electrostatic potentials the need for high sampling density in formamide conformational analysis. J. Comp. Chem. 1990, 11, 361-373. [26] Singh, U. C.; Kollman, P. A. An approach to computing electrostatic charges for molecules. J. Comp. Chem. 1984, 5, 129-45. [27] Besler, B. H.; Merz Jr, K. M.; Kollman, P. A. Atomic charges derived from semiempirical methods. J. Comput. Chem. 1990, 11, 431-439. [28] Verstraelen, T.; Ayers, P. W.; van Speybroeck, V.; Waroquier, M. Hirshfeld-E Partitioning: AIM Charges with an Improved Trade-off between Robustness and Accurate Electrostatics. J. Chem. Theory Comput. 2013, 9, 2221-2225. [29] Bader, R. W. F.; Nguyen-Dang, T. T. Quantum Theory of Atoms in MoleculesDalton Revisited. Adv. Quantum Chem. 1981, 14, 63-124. [30] Hirshfeld, F. L. Bonded-atom fragments for describing molecular charge densities. Theoret. Chim. Acta. 1977, 44, 129-138. [31] Storer, J.; Giesen, D.; Cramer, C.; Truhlar, D. G. Class IV charge models: a new semiempirical approach in quantum chemistry. J. Comput.-Aided Mol. Des. 1995, 9, 87-110. [32] Li, J.; Zhu, T.; Cramer, C. J.; Truhlar, D. G. New Class IV Charge Model for Extracting Accurate Partial Charges from Wave Functions. J. Phys. Chem. A. 1998, 102, 1820-1831. [33] Thompson, J. D.; Cramer, C. J.; Truhlar, D. G. Parameterization of charge model 3 for AM1, PM3, BLYP, and B3LYP. J. Comput. Chem. 2003, 24, 12911304. [34] Kelly, C. P.; Cramer, C. J.; Truhlar, D. G. SM6: A Density Functional Theory Continuum Solvation Model for Calculating Aqueous Solvation Free Energies of Neutrals, Ions, and Solute-Water Clusters. J. Chem. Theory Comput. 2005, 1, 1133-1152. ACS Paragon Plus Environment 22

Page 22 of 35

Page 23 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


[35] Olson, R. M.; Marenich, A. V.; Cramer, C. J.; Truhlar, D. G. Charge Model 4 and Intramolecular Charge Polarization J. Chem. Theory Comput. 2007, 3, 2046-2054. [36] Marenich, A. V.; Jerome, S. V.; Cramer, C. J.; Truhlar, D. G. Charge Model 5: An Extension of Hirshfeld Population Analysis for the Accurate Description of Molecular Interactions in Gaseous and Condensed Phases. J. Chem. Theory Comput. 2012, 8, 527-541. [37] Dodda, L. S.; Vilseck, J. Z.; Cutrona, K. J.; Jorgensen, W. L. Evaluation of CM5 Charges for Nonaqueous Condensed-Phase Modeling. J. Chem. Theory Comput. 2015, 11, 4273-4282. [38] Cole, D. J.; Vilseck, J. Z.; Tirado-Rives, J.; Payne, M. C.; Jorgensen, W. L. Biomolecular force field parameterization via Atoms-in-Molecules Electron Density Patitioning. J. Chem. Theory Comput. 2016, 12, 2312-2323. [39] Macchiagodena, M.; Mancini, G.; Pagliai, M.; Barone, V. Accurate prediction of bulk properties in hydrogen bonded liquids: amides as case studies. Phys. Chem. Chem. Phys. 2016, 18, 25342-25354. [40] Jorgensen, W. L.; Tirado-Rives, J. The OPLS [optimized potentials for liquid simulations] potential functions for proteins, energy minimizations for crystals of cyclic peptides and crambin. J. Am. Chem. Soc. 1988, 110, 1657-1666. [41] Virtual Chemistry. http://www.virtualchemistry.org/ (accessed November 15, 2012). [42] Gaussian 09, Revision B.01, 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. Petersson, H. Nakatsuji, M. Caricato, X. Li, H. P. Hratchian, A. F. Izmaylov, J. Bloino, G. Zheng, J. L. Sonnenberg, M. Hada, M. Ehara, 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, F. Ogliaro, M. Bearpark, J. J. Heyd, E. Brothers, K. N. Kudin, V. N. Staroverov, T. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. Rendell, 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, 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, G. A. Voth, P. Salvador, J. J. Dannenberg, S. Dapprich, A. D. Daniels, O. Farkas, J. B. Foresman, J. V. Ortiz, J. Cioslowski, and D. J. Fox, Gaussian, Inc., Wallingford CT, 2010.



[43] Zhao, Y.; Truhlar D. G. The M06 suite of density functionals for main group thermochemistry, thermochemical kinetics, noncovalent interactions, excited states, and transition elements: two new functionals and systematic testing of four M06-class functionals and 12 other functionals. Theor. Chem. Acc. 2008, 120, 215-241. [44] McLean, A. D.; Chandler, G. S. Contracted Gaussian-basis sets for molecular calculations. 1. 2nd row atoms, Z=11-18. J. Chem. Phys. 1980, 72, 5639-5648. [45] Krishnan, R.; Binkley, J. S.; Seeger, R.; Pople, J. A. Self-Consistent Molecular Orbital Methods. 20. Basis set for correlated wave-functions. J. Chem. Phys. 1980, 72, 650-654. [46] Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for highly efficient, load-balanced, and scalable molecular simulation. J. Chem. Theory Comput. 2008, 4, 435-447. [47] Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. G. E. M. LINCS: A linear constraint solver for molecular simulations. J. Comput. Chem. 1997, 18, 14631472. [48] Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen L. G. A smooth particle mesh Ewald method. J. Chem. Phys. 1995, 103, 8577-8593. [49] Hansen, J. P.; McDonald, I, R. Theory of Simple Liquids. 2013, 4th edition, Elsevier Ltd. [50] Neumann, M. Dipole moment fluctuation formulas in computer simulations of polar systems. Mol. Phys. 1983, 50, 841-845. [51] Truckymchuk, A.; Alejandre, J. Computer simulations of liquid/vapor interface in Lennard-Jones fluids: Some questions and answers J. Chem. Phys. 1999, 111, 8510-8523. [52] Orea, P.; López-Lemus, J.; Alejandre J. Oscillatory surface tension due to finitesize effects. J. Chem. Phys. 2005, 123, 114702-114813. [53] Gonzalez-Melchor, M.; Bresme, F.; Alejandre, J. Molecular dynamics simulations of the surface tension of ionic liquids. J. Chem. Phys. 2005, 122, 104710. [54] Lide, D.R.; Haynes, W.M. CRC Handbook of Chemistry and Physics, 90th ed. CRC Press, Boca Raton, FL, 2009. [55] Payne, R.; Theodorou, I. E. Dielectric Properties and Relaxation in Ethylene Carbonate and Propylene Carbonate. J. Phys. Chem. 1972, 76, 2892-2900. ACS Paragon Plus Environment 24

Page 24 of 35

Page 25 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


[56] Xie, W.; Pu, J.; Mackerell, A. D.; Gao. J. Development of a Polarizable Intermolecular Potential Function (PIPF) for Liquid Amides and Alkanes. J. Chem. Theory Comput. 2007, 3, 1878-1889. [57] Yaws C. L. Thermophysical properties of chemicals and hydrocarbons, William Andrew, 1999. [58] Wohlfarth, Ch.; Wohlfarth, B. Surface Tension of Pure Liquids and Binary Liquid Mixtures, Simmetry Springer, New York, NY, 1994. [59] Beesley, A. H.; Evans, D. F.; Laughlin, R. G. Evidence for essential role of hydrogen bonding in promoting amphiphilic self-assembly: Measurement in 3methylsydnome. J. Phys. Chem. 1988, 92, 791-793. [60] Nayak, J. N.; Aralaguppi, M. I.; Toti, U. S.; Aminabhavi, T. M. Density, Viscosity, Refractive Index, and Speed of Sound in the Binary Mixtures of Trin-butylamine + Triethylamine, + Tetrahydrofuran, + Tetradecane, + Tetrachloroethylene, + Pyridine, or +Trichloroethylene at (298.15, 303.15, and 308.15) K. J. Chem. Eng. Data. 2003, 48, 1483-1488. [61] Enders, S.; Kahl, H.; Winkelmann, J. Surface tension of the ternary system water + acetone + toluene. J. Chem. Eng. Data. 2007, 52, 1072-1079. [62] Barthel, J.; Neueder, R.; Roch, H. Density, Relative Permittivity, and Viscosity of Propylene Carbonate + Dimethoxyethane Mixtures from 25 ◦ C to 125 ◦ C. J. Chem. Eng. Data. 2000, 45, 1007-1011. [63] LaPlanche, L. A.; Rogers, M. T. cis and trans Configurations of the Peptide Bond in N-Monosubstituted Amides by Nuclear Magnetic Resonance. J. Am. Chem. Soc. 1964, 86, 337-341. [64] Verevkin, S. P.; Toktonov, A. V.; Chernyak, Y.; Sch¨ affner, B.; B¨ orner, A. Vapour pressure and enthalpy of vaporization of cyclic alkylene carbonates. Fluid Phase Equilibria. 2008, 268, 1-6. [65] Norinaga, K.; Wargardalam, V. J.; Takasugi, S.; Iino M.; Matsukawa, S. Measurement of Self-Diffusion Coefficient of Asphaltene in Pyridine by Pulsed Field Gradient Spin-Echo H NMR. Energy Fuels. 2001, 15, 1317-1318. [66] Guevara-Carrion, G.; Janzen, T.; Mu˜ noz-Mu˜ noz Y. M.; Vrabec, Y. Mutual diffusion of binary liquids mixtures containing methanol, ethanol, acetone, benzene, cyclohexane, toluene and carbon tetrachloride. J. Chem. Phys. 2016, 144, 124501.



[67] Das, K. P.; Ceglie, A.; Lindman, B. Microstructure of formamide microemulsions from NMR self-diffusion measurements. J. Phys. Chem. 1987, 91, 2938-2946. [68] Della-Volpe, C.; Guarino, G; Sartorio, R.; Vitaliano, V. Diffusion, Viscosity, and Refractivity Data on the System Dimethylformamide-Water at 20 and 40 o C. J. Chem. Eng. Data. 1986, 31, 37-40. [69] Kawanishi, H.; Tsunashima, Y.; Horii F. A nest of structures in dynamics of cellulose diacetate in N,N-dimethylacetamide in quiescent solution state studied by dynamic light scattering. J. Chem. Phys. 1998, 109, 11027. [70] O’Reilly, D. E.; Peterson E. M. Self-Diffusion Coefficients and Rotational Correlation Times in Polar Liquids. II. J. Chem. Phys. 1971, 55, 2155. [71] Elliot, A. J.; McCracken D. R.; Buxton, G. V.; Wood, N. D. Estimation of Rate Constants for Near-diffusion-controlled Reactions in Water at High Temperatures. J. Chem. Soc. Faraday Trans. 1990, 86, 1539-1547. [72] Holz, M.; Mao, X.; Seiferling, D. Experimental study of dynamic isotope effects in molecular liquids: Detection of translation-rotation coupling. J. Chem. Phys. 1996, 104, 669. [73] Kondo, K.; Sano, M.; Hiwara, A.; Omi, T.; Fujita, M.; Kuwae, A.; Iida, M.; Mogi, K.; Yokoyama, H. Conductivity and Solvation of Li+ Ions of LiPF6 in Propylene Carbonate Solutions. J. Phys. Chem. B. 2000, 104, 5040-5044. [74] Cases, A. M.; Gómez Margliano, A. C.; Bonatti, C. M.; Sólimo, H. N. Density, Viscosity, and Refractive Index of Formamide, Three Carboxylic Acids, and Formamide + Carboxylic Acid Binary Mixtures. J. Chem. Eng. Data. 2001, 46, 712-715. [75] Nikam, P. S.; Jadhav, M. C.; Hasan, M. Density and Viscosity of Mixtures of Nitrobenzene with Methanol, Ethanol, Propan-1-ol, Propan-2-ol, Butan-1-ol, 2Methylpropan-1-ol, and 2-Methylpropan-2-0l at 298.15 and 303.15 K. J. Chem. Eng. Data. 1995, 40, 931-934. [76] Chein-Hsiun, T.; Shu-Lien, L.; I-Hung, P. Excess Volumes and Viscosities of Binary Mixtures of Aliphatic Alcohols (C1-C4) with Nitromethane. J. Chem. Eng. Data. 2001, 46, 151-155. [77] Pal, A.; Kumar, A. Excess Molar Volumes, Viscosities, and Refractive Indices of Diethylene Glycol Dimethyl Ether with Dimethyl Carbonate, Diethyl Carbonate, and Propylene Carbonate at (298.15, 308.15, and 318.15) K. J. Chem. Eng. Data. 1998, 43, 143-147. ACS Paragon Plus Environment 26

Page 26 of 35

Page 27 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


[78] Nayak, J. N.; Aralaguppi, M. I.; Toti, U. S.; Aminabhavi, T. M. Density, Viscosity, Refractive Index, and Speed of Sound in the Binary Mixtures of Trin-butylamine + Triethylamine, + Tetrahydrofuran, + Tetradecane, + Tetrachloroethylene, + Pyridine, or + Trichloroethylene at (298.15, 303.15, and 308.15) K. J. Chem. Eng. Data. 2003, 48, 1483-1488. [79] Tsierkezos, N. G. Cyclic Voltammetric Studies of Ferrocene in Nonaqueous Solvents in the Temperature Range from 248.15 to 298.15 K. J Solution Chem. 2007, 36, 289-302. [80] Manzanilla-Granados, H. M. ; Saint-Mart´ın, H.; Fuentes-Azcatl, R.; Alejandre, J. Direct Coexistence Methods to Determine the Solubility of Salts in Water from Numerical Simulations. Test Case NaCl. J. Phys. Chem. B. 2015, 119, 8389-8396. [81] Jelinska-Kazimierczuk, M.; Szydlowski, J. Physicochemical Properties of Solutions of Amides in H2 O and in D2 O. J. Solution Chem. 2001, 30, 623-640. [82] Papamatthaiakis, D.; Aroni, F.; Havredaki, V. Isentropic compressibilities of (amide + water) mixtures: A comparative study. J. Chem. Thermodyn. 2008, 40, 107-118. [83] Rohdewald, P.; Moldner, M. Dielectric Constants of Amide-Water Systems. J. Phys. Chem. 1973, 77, 373-377. [84] Sengwa, R. J.; Khatri, V.; Sankhla, S. Static Dielectric Constants and Kirkwood Correlation Factor of the Binary Mixtures of N-Methylformamide with Formamide, N,N-Dimethylformamide and N,N-Dimethylacetamide. J. Solution Chem. 2009, 38, 763-769. [85] Sengwa, R. J.; Khatri, V.; Sankhla, S. Static Dielectric Constant, Excess Dielectric Properties, and Kirkwood Correlation Factor of Water-Amides and WaterAmines Binary Mixtures. Proc. Indian Natn. Sci. Acad. 2008, 74, 67-71. [86] Wensink, E. J. W.; Hoffmann, A. C.; van Maaren, P. J.; van der Spoel, D. Dynamic properties of water/alcohol mixtures studied by computer simulation. J. Chem. Phys. 2003, 119, 7308-7317. [87] Zaichikov, A. M. Excess Enthalpies of Binary Mixtures of Formamide with Secondary Amides. J. Thermal Anal. 1998, 54, 279-284. [88] Zielkiewicz, J. Excess volumes in (N-methylformamide+methanol+water) at the temperature 313.15 K. J. Chem. Thermodyn. 1995, 27, 1275-1279.


Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58

ACS Paragon Plus Environment

Page 28 of 35

Page 29 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41




g(r)(NH...OH)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38

8

Page 30 of 35

Acetamide Formamide N-Methylformamide

6

4

2

0 0.1

0.2


0.3 r[nm]

0.4

0.5

Page 31 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41




10 8 6 4 2 0 -2

200

Number of clusters

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

150

0

5

10

15

Page 32 of 35

20

25

100 Acetamide Formamide N-methylformamide 50

0

0

5

10 15 Zize of clusters


20

25

Page 33 of 35

1300

1060

1075 B)

C)

A) 1050

Exp NewFF OPLS/AA GAFF CGenFF

1250

1040

1050

1200

3

ρ [kg/m ]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47


1030 1150

1025

1020 1100 1010 1000 1050 1000

990

1000 0

10 20 X(ACE)

30

0

20

40 60 X(FM)


80

100

975

0

20

40 60 80 X(NMFM)

100

100

125


A) 90

B)

100

80 75

ε

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36

Page 34 of 35

70 50

60 50

0

10

250 C)

200

30

40

Exp Exp NewFF OPLS/AA GAFF CGenFF

25

0

200

20

40 60 X(FM)

80

100

40 60 X(NMF)

80

100

D)

150

ε

150

20 X(ACE)

100

100

50

50 0

0

20

40 60 X(NMF)

0


80

100

0

20

3

Page 35 of 35

B)

A) 0 2

-1

3

∆VMix.[cm /mol]

∆HMix.[kJ/mol]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44


1

-2

0

Exp. NEWFFP OPLS/AA GAFF CGenFF

-3

-1

0

0.2

0.4 0.6 X(NMF)

0.8

-4


1

0

0.2

0.4 0.6 X(NMF)

0.8

1

Force field parameterization from the Hirshfeld molecular electronic

Recommend Documents