Structure and Weak Hydrogen Bonds in Liquid Acetaldehyde

Investigaram-se a estrutura e formação de pontes de hidrogênio no acetaldeído líquido por simulações de Monte Carlo. Um modelo para acetaldeído com todos os átomos explicitados foi otimizado no presente trabalho. Os valores teóricos obtidos para a entalpia de vaporização e a densidade do líquido mostram boa concordância com os dados experimentais. Os gráficos de função de distribuição radial (rdf) indicam um líquido bem estruturado comparado com outros líquidos orgânicos dipolares semelhantes. A minimização usando mecânica molecular em fase gasosa produz um trímero com uma estrutura muito estável. A geometria deste complexo apresenta boa concordância com os gráficos de rdf. A correlação sítio-sítio mais próxima observada é entre o oxigênio e o hidrogênio da carbonila, sugerindo que esta correlação possui um importante papel na estrutura do líquido e suas propriedades. A distância média O⋅⋅⋅H e o ângulo C–H⋅⋅⋅O obtidos são característicos de pontes de hidrogênio fracas.


Introduction
2][3] The development of suitable computer simulation methods have made it feasible to investigate a great number of liquid state properties that would not be possible, or difficult to access, with other theoretic or experimental techniques. 4,57][8][9] Therefore, the optimisation of meaningful intermolecular potential functions is of fundamental importance for the investigation of liquid state properties using theoretical methods.This optimisation is not simple, however, due to the complexity of the liquid intermolecular interactions, and some criteria, such as pairwise additivity and consideration of dipolar interaction exclusively, have been used to simplify the interaction potential making the work feasible. 10,11e have recently contributed to the development of parameterised potential functions to study dipolar organic liquids, using these approximations. 12,133][14][15] The study has been centred in amides, since these molecules are suitable models for peptide moieties, and dimethyl sulfoxide because these liquids are popular solvents in organic chemistry with stressed effects in reaction kinetics. 16In this work, new additive pairwise potential has been parameterised for acetaldehyde.These results are part of a more extensive work aiming to investigate hydrogen bonding in biologically relevant systems, and extend the data earlier published. 12,14In a future paper, thermodynamics results for aqueous solutions of acetaldehyde, formamide and methylformamide will be reported.0][21] Gas phase molecular mechanics (MM) have been used to improve the insights about the species prevailing in the liquid phase.

Methodology
The intermolecular potential function Following usual procedures in force field calculations, 22 the molecules were modelled by collections of interacting sites.The energy E ab between molecules a and b was given by a sum of electrostatic Coulomb plus nonbonded Lennard-Jones terms centred on the sites: (1)   where r ij is the distance between site i in a and site j in b and q i and q j are fractional point charges located on the i and j molecular sites.For each site k, the parameters A kk and B kk were given by A kk =4ε k σ k 12 and B kk =4ε k σ k 6 , where ε k and σ k are the Lennard-Jones parameters for the k-th site.In most of the systems (mainly in the case of polarizable species), many-body terms can play an important role in the interaction potential.However, despite the weight of these terms in the potential, they are rarely included in the computer simulations, for this would be very expensive.The average many-body effects can be partially included in the pairwise approximation by defining an effective pair potential.
Parameters A ij and B ij for a nondiagonal interaction [i,j] were obtained using the geometric combining rules 23 A ij = (A ii A jj ) 1/2 and B ij = (B ii B jj ) 1/2 .This is a standard procedure used in this field of study.In an earlier work it was shown that results obtained using the combining rule above and the equally well-known Lorentz-Berthelot rule are equivalent. 24

Molecular model and potential parameters
The acetaldehyde optimised in this work is a rigid allatom model, the geometry of which was taken from allatom OPLS model. 25Point charges needed to take into account the Coulomb interaction in equation (1) were obtained from Mulliken charges calculated at MP2 level with a 6-31G* basis set.The calculations were performed with the Gaussian 94 program. 26The utilisation of CHELPG procedure 27 at the RHF level to calculate the point charges has been suggested in an earlier work. 28Extensive studies have revealed, however, that for molecules with a carbonyl group, mainly those presenting peptide bond, Mulliken charges better represent the molecular charge asymmetry than the CHELPG procedure. 29Mulliken population analysis was used to parameterise the atom partial charges in CHARMM program. 9It has also been suggested that charges obtained using CHELPG are not dependent on the inclusion of electron correlation. 30However, previous calculations revealed that RHF and MP2 charges are substantially different and electron correlation must be included in the charges calculations of this type of molecules. 29he Lennard-Jones potential parameters were optimised using the all-atom OPLS parameters as a starting trial set. 25They were then adjusted through MC interactive simulations to provide a satisfactory fitting for experimental results of density and enthalpy of vaporisation of acetaldehyde at 18 °C and 1 atm.This adjustment was performed manually.This procedure is preferable over automatic ones, since they often lead to a combination of unphysical parameters that reproduce the input data. 9

Monte Carlo simulations
Thermodynamic and structural properties have been computed for the pure liquid by MC simulations. 31etropolis importance sampling and periodic boundary conditions 17 have been used in the NPT ensemble.Cubic boxes containing 500 molecules were simulated at 1 atm and 18 °C.Starting from a molecular distribution in a face centred cell, new configurations were generated by randomly rotating and translating a molecule along randomly chosen Cartesian directions.The ranges for monomer translations and rotations were established to provide an acceptance/trial ratio of about 45% for new configurations.A volume movement was attempted at every 500 configurations and the range for volume change was chosen to provide the same acceptance/trial ratio as above.After a volume movement the intermolecular distances were scaled suitably.
The computational cost of calculating all nonbonded interactions in a system scales as the square of the number of atoms and becomes expensive for a large number of them.It is common practice to use a cut off beyond which these interactions are neglected.In this work, to calculate the intermolecular energy (equation 1), a full intermolecular interaction was considered whenever any site-to-site distance was lower than a cut-off radius r c = 11.0Å.The contributions of Lennard-Jones interactions beyond the cut-off radius were taken into account using the methodology presented elsewhere. 23In this methodology the long-range contribution from the Lennard-Jones potential is calculated by integrating the potential function from the cut-off radius to infinity, assuming that the pair correlation function is unitary.On the other hand, as it has been widely used in this field of study, long-range corrections to take into account Coulombic interactions beyond cut-off radius were not included in the present work.A brief analysis on the aspects that validate this assumption was done recently. 32When the systems under investigation are constituted of neutral molecules, the main long-range contribution from the electrostatic potential to the cohesive energy comes from the dipole-dipole interaction.This contribution decrease asymptotically with 1/r 3 , where r is roughly the distance between the centre of mass of two interacting molecules and therefore the dipoledipole interaction averaged over the angular orientation decrease asymptotically with 1/r 6 . 33Consequently, the long-range electrostatic contribution to the electrostatic energy is lower than could be expected for a strictly Coulombic system.This is in accord with recent results 34 using the generalised Born model. 35,36Also in good agreement with this, it is observed that site-site pair correlation functions, thermodynamic properties and single particle dynamics of polar liquids, appears to be quite insensitive to the long-range forces in the system. 21he calculations started with an equilibration phase with 1.6x10 7 configurations and the averages were then obtained from a new segment with other 1.6x10 7 configurations.Statistical uncertainties were calculated from separate averages over blocks of 2x10 5 configurations.

Potential parameters for acetaldehyde
Table 1 lists the final set of potential parameters obtained for acetaldehyde.The parameterisation of a new set of parameters for acetaldehyde is justified because the OPLS charge for aldehydic hydrogen (zero), does not fit properly the weak hydrogen bond features presented for molecules like acetaldehyde. 37The study of these weak hydrogen bonds is the main aim of our work.
It is usual in liquid phase simulations in which molecular dynamic polarisation effects have not been included, to consider as desirable a theoretical molecular dipole moment around 15-20 % greater than the experimental one in the gas phase. 28,38,39This overestimation in the dipole moment has been used as a manner to compensate for neglecting polarisation effects in the simulations.However, the dipole moment for the acetaldehyde molecular charge distribution optimised in this work is µ = 2.6 D, very similar, therefore, to the experimental value µ = 2.7 D. 40 Theoretical values of thermodynamic properties for acetaldehyde obtained with OPLS and the present model are shown in Table 2.The experimental values are also reported for comparison.

Radial correlation function and the liquid structure
Radial distribution functions (rdf) play an important role in the understanding of the liquid structure.Experimental information about rdf can be obtained by neutron scattering and X-ray diffraction.Nevertheless, the large number of overlapping distances, both within and between molecules, and the large number of different sitesite correlation functions added together in the diffraction pattern, make the experimental investigation of liquid structure difficult. 32Therefore, the use of molecular modeling and computer simulations to acquire the various site-site correlation functions, can be useful in providing improved insight about the liquid structure.
Figure 1 shows the rdf computed with the present model for acetaldehyde.The broadness and positions of the peaks indicate that acetaldehyde is a liquid relatively well structured compared to other organic liquids like tetrahydrofuran and acetone, 24,32 which present site-site correlation peaks above 3.0 Å and broad peaks instead of the quite sharp peaks observed in the present case.2][43][44] The distance O/H(C=) is characteristic of weak hydrogen bonds. 45The same can be addressed for the correlation between oxygen and the methyl hydrogens.
As it can be noticed, 2 hydrogens of the methyl group are identically correlated with oxygen while the third is further apart.The shortest O/H distance is about 2.6 Å which permit to classify this interaction also as a weak hydrogen bond, using the same criterion used for the earlier case.The features of weak hydrogen bonds and the criteria used to classify them have been extensively discussed very recently by many authors [46][47][48][49][50][51][52][53][54][55][56][57] and have already been used by us in previous works. 12,14It is important to recall that hydrogen bonds are additive interactions 45 and, because of this characteristic, even weak hydrogen bonds can play an important role in stabilising the structure of a material.This can be the reason behind the observation that the methyl carbon is closer to oxygen than the carbonyl carbon.This aspect, however, deserves a closer attention, since other aspects are involved in these interactions.One has to take into account that the methyl carbon is 0.25 Å smaller in diameter than the carbonyl one, and the lower distance verified between oxygen and methyl carbon compared to the distance between oxygen and the carbonyl carbon could be just consequence of this fact.On the other hand, the charge on the methyl carbon is quite negative, as the charge of the oxygen, while the charge of the carbonyl carbon is positive.Considering Coulombic effects, it would be expected an approximation of oxygen to the carbonyl carbon, due to the attraction between them and the repulsion between oxygen and the methyl carbon.However, the contrary behaviour is observed.In this scenario, one could suppose that the attractive interaction between the methyl hydrogens and the oxygen is responsible for this found, but more details to this respect are not possible to be addressed from the present simulations.In this sense, higher level calculations could be interesting to clarify this issue.
The O/C(me) and O/C= pair correlations present some interesting features that is worth being discussed.The O/C= correlation function presents two peaks, one of them is in a position almost overlaping the peak for O/C(me) correlation.This means that the oxygen is roughly localised in the bissecting line of the angle OCC of neighbouring molecules.The second peak about 4.4 Å as well as some shoulders in the O/H(C=) and O/O rdf suggest that the bulk of the liquid should contain more complicated species than just "dimers".Molecular mechanics at MM2 level 58 was used in order to minimise a structure that could explain the rdf's shown in Figure 1.After the corresponding calculations it was possible to rise the "trimer" shown 61 in Figure 2. The molecules are not coplanar and there is an extraordinary agreement between the geometry of the trimer and the plots of g(r) shown in Figure 1.Little differences in the distances between the same sites taking pairs of molecules (tenth of Å) are found, which explains the pattern observed in the pair correlation peaks.However, the main feature associated with that trimer is that there are two O/C= pairs located at distances that closely agree with the O/C= pair rdf shown in the Figure 1.A standard procedure in this field is the calculation of the site-site coordination number by properly integrating the correspondent correlation function over the first maximum up the position of the first minimum. 23The coordination number N r 0 is given for , where ρ is the liquid density and r is the integrating radius.r 0 is the radius at the first minimum.
The number of hydrogen bonds O/H(C=) obtained was 1.8.This result corroborates that a trimer is the representative structure in the liquid bulk.Thus, each one of the oxygens in the trimer complex is near to two equidistant carbonilic hydrogens.The C/H(C=) distance calculated in the complex also compares well with the corresponding site-site distance in the rdf graphic.][43][44][45][46][47][48][49][50][51][52][53][54][55][56][57] Besides, the relative position of one molecule with respect to an other, due to the correlation between the oxygen and the methyl hydrogens, could be interpreted as more a indication of the weight of the weak hydrogen bonds in the structure.Taking this result into account one can suppose that this type of interaction can play an important role in the structure of molecules like proteins for example, mainly if normal hydrogen bonds are not involved.
The energy of the acetaldehyde trimer was also investigated using molecular mechanics at MM2 level, 58 as it was done in a previous work for acetone. 32The acetaldehyde trimer here reported has a stabilisation energy at MM2 level of -10.23 kcal/mol, whilst the stabilisation energy of the acetone dimer 32 is -2.91 kcal/mol.Thus the trimer of acetaldehyde seems to be particularly more stable than the dimer of acetone.The O/H interactions can be playing a considerable role in this stabilisation.The high energy of stabilisation of the acetaldehyde trimer permits to understand why the liquid is so well structured compared to acetone, according to the rdf's of both of them.
Finally, the rdf obtained with the acetaldehyde molecule parameterised in this work are quite similar to those calculated with the OPLS model despite the huge difference in the atomic charges of the two models.The similarity of their dipole moment, however, would indicate that the dipole moment would play an important role in the liquid structure, as have already been pointed out previously. 59According to Amey, 59 there is a preferential inter-molecular geometry arising from dipole-dipole orientation, and any form of molecular array in the liquid must be related to the dipole moment.On the other hand, Luzar et al. 60 have shown that a substantial part of the sitesite correlation in DMSO is due to packing effects and, hence, cannot be entirely credited to molecular association.Taking into account these results it is difficult to quantify how much the liquid is structured because of the dipoledipole/packing effects and how much because of atomic interactions.A theoretical investigation aiming to quantify these effects would be extremely important.Of course, even if the molecules were relatively oriented because of the dipolar interaction, specific atomic roles would complete the scenario giving the final structure of the liquid.

Conclusions
A new model for acetaldehyde has been optimised, and some thermodynamic and structural features of it have been calculated using Monte Carlo simulations.The site charges used, obtained from Mulliken distribution at MP2 level, seem to provide a good charge distribution for the molecule.The charge distribution fits well the experimental value the dipole moment.The model here reported has a suitable ability to mimic the energetic and structural properties of liquid acetaldehyde that are related to density and heat of vaporisation.The rdf for site-site correlations present quite sharp peaks indicating that acetaldehyde is a well organised liquid.Molecular mechanics calculations showed that the predominant species in the bulk is a trimer, whose structure was here reported.This trimer seems to be quite stable.The relative positions among atoms in this complex suggest the formation of O/carbonyl hydrogen and O/methyl hydrogens weak hydrogen bonds.These interactions play an important role in the stabilisation of the liquid structure.

Table 1 .
Potential parameters for acetaldehyde