Article pubs.acs.org/JPCB

Molecular Dynamics Simulation Study of Methanesulfonic Acid Manel Canales*,† and Carlos Alemán‡ †

Departament de Física i Enginyeria Nuclear, Universitat Politècnica de Catalunya, Campus Nord-Edifici B4-B5, Jordi Girona 1-3, Barcelona E-08034, Spain ‡ Departament d’Enginyeria Química, Universitat Politècnica de Catalunya, Diagonal 647, Barcelona E-08028, Spain ABSTRACT: A molecular dynamics simulation study of methanesulfonic acid has been carried out using a reliable force field in a large range of temperatures. Several thermodynamic, structural, and dynamical properties have been calculated and compared with the available experimental data. The density, the shear viscosity, the heat of vaporization, and the melting temperature results, calculated from this force field, are in a good agreement with the experimental data. Analysis of the influence of the hydrogen bonds in structural and dynamical properties has also been performed. The continuous and interrupted methodologies to compute hydrogen bonding lifetimes have been applied. The interrupted hydrogen bond lifetimes values are consistent with the diffusion and viscosity coefficients. The activation energies of the self-diffusion, the reorientational motions, and the hydrogen bonding lifetimes are coincident.



INTRODUCTION Methanesulfonic acid (CH3SO2OH, MSA), which is liquid at room temperature, is the simplest of the sulfonic acids.1 Several physical and chemical properties, such as its nonoxidizing nature, strong acidity, high solubility of its salts, low toxicity, thermal stability, and biodegradability, are the keys to explain the wide industrial use of MSA as solvent and reagent.2 For example, it is used in electroplating,2 in acid catalysis for esterification, polymerization, and alkylation reactions,3 in electronics to increase the conductivity of conjugated polymers,4 etc. Moreover, it plays an important role in biological and environmental systems. Thus, MSA is an atmospheric oxidation product of the dimethyl sulfide produced by marine biota, and therefore, it is an indicator of the atmospheric sulfur cycle.5,6 Computer simulation is a useful tool that can be employed to validate models of physical systems by comparing the simulation results with available experimental data. The quality of a model will be established by its ability to predict the maximum number of properties. In the particular case of MSA, several ab initio studies, using the Gaussian package, have been performed to compare enthalpies of formation, structures, and spectroscopic properties determined experimentally and theoretically.7−9 Similar studies have also been made to study the structure and energy of hydrated MSA clusters10−12 and the equilibrium structure of MSA-pyridine.13 However, to our knowledge, a very scarce number of MSA molecular dynamics (MD) studies can be found in the literature.14 Therefore, the first objectives of this work are, on one hand, to investigate MSA at the atomistic level in a large range of temperatures and, on the other hand, to develop a reliable force field by © 2014 American Chemical Society

comparing the simulations results with the available experimental data. The understanding of the hydrogen bond interaction is of paramount importance to explain many of the properties of a wide variety of chemical and biological systems.15,16 In this regard, it is decisive to know details such as, for instance, the mean number of hydrogen bonds per molecule, the lifetime, and its strength. Computer simulations of water and other hydrogen bonding liquids, such as alcohols, have been performed during the past four decades.17 As it is wellknown, liquid water18 is characterized by a topologically complex three-dimensional hydrogen-bonded arrangement, which, on average, resembles the tetrahedral network of hydrogen bonds that exists in the case of ice. However, for methanol there is only one donor hydrogen and a methyl group that force the system to organize in a series of nonlinear hydrogen-bonded open chains without significant branching.19 MSA is characterized by the presence of two types of oxygen atoms: the sulfonyl ones, which are bonded to the sulfur by a double bond, and the hydroxyl oxygen atoms, which are bonded to the sulfur and the hydrogen by single bonds. Thus, the ultimate goal of this work is to study the influence of hydrogen bonds in structural and dynamical properties analyzing the role of both types of oxygen atoms. The paper is organized as follows. Simulation details are given in the next section. Thermodynamic, structural and dynamical properties, and hydrogen bond analyses are Received: January 23, 2014 Revised: March 3, 2014 Published: March 4, 2014 3423

dx.doi.org/10.1021/jp500817s | J. Phys. Chem. B 2014, 118, 3423−3430

The Journal of Physical Chemistry B

Article

have been computed using the set of partial charges collected in Table 4. These values have been employed by other authors to study similar systems. Thus, the charges associated with the atoms of the CH3 group are those of the OPLS-AA (all atom) model of Jorgensen et al.26 and have been used by Taylor and Shields27 to study the ethanol liquid−vapor interface. The charges of the SO2OH group are those employed by Ding et al.28 to analyze the distribution of sulfuric acid clusters and by Paddison29 to study trifluoromethanesulfonic acid−water clusters. MD simulations in the isothermal−isobaric (NPT) ensemble of a system with 500 MSA molecules, using the usual periodic boundary conditions, with an integration time step of 1 fs have been performed. The van der Waals interactions have been truncated at the cutoff distance of 1 nm. Electrostatic interactions were computed by means of Ewald summations.30 The real space term was defined by the van der Waals cutoff, while the reciprocal space was computed by interpolation into an infinite grid of points (particle mesh Ewald) with maximum space grid points being 0.12 nm.31 Before run the production MD trajectory, the starting structure was equilibrated using the following strategy: a first configuration has been generated introducing 500 MSA molecules at random inside a cubic box with sides of 4 nm. Next, 104 steps of energy minimization, using a steepest descent algorithm,32 have been applied to relax the system. Finally, 10 ns of NPT-MD at 298 K and 1 atm was used to equilibrate the system. Both temperature and pressure were controlled by the weak coupling method of the Berendsen thermobarostat,33 using a time constant for the heat bath coupling of 0.1 ps and a pressure relaxation time of 1 ps. The dynamical and structural properties have been calculated from the configurations generated during 107 NPT-MD time steps, with the coordinates and velocities being saved every 100 steps. Finally, to study the temperature dependence of some of these properties, a set of 21 analogous NPT-MD simulations (at 1 atm) in the range between 200 and 400 K with a temperature step of 10 K have also been performed.

presented in the Results and Discussion section. Finally, some concluding remarks are gathered in the last section.



METHODOLOGY The Gaussian 03 program20 has been used to determine the equilibrium geometry of a single MSA molecule. In particular, the bond length, the bending bond angle, and the torsional dihedral angle have been computed. Calculations have been performed at the density functional level of theory (DFT) using the Becke’s three-parameter exchange hybrid functional (B3)21 and the nonlocal correlation functional of Lee, Yang, and Parr (LYP).22 This exchange-correlation functional (B3LYP) was combined with the 6-31+G(d,p) basis set. MD simulations have been carried out using the GROMACS 4.6 software package. 23 The bonded and nonbonded interactions have been calculated according to the GROMOS force field functional form.24 The bond stretching distance, bond angle bending, and torsional dihedral angle parameters have been computed using, on one hand, bond lengths, bond angles, and torsional angles calculated from our DFT calculations and, on the other hand, the force constant parameters of the GROMOS 53 A6 force field version.25 All these parameters are gathered in Tables 1, 2, and 3. The Table 1. Bond Lengths and Stretching Constants bond

length (nm)

Kb (106 kJ mol−1 nm−4)

S−C S−O′ SO

0.1795 0.1653 0.1465

5.94 8.37 8.37

bond

length (nm)

Kb (106 kJ mol−1 nm−4)

O′−H′ C−H

0.0973 0.1090

15.7 12.3

Table 2. Bond Angles and Harmonic Force Constants bond

angle (deg)

Kθ (kJ mol−1)

H′−O′−S O′−S−O O−S−O O′−S−C

110 110 120 100

550 450 450 503

bond

angle (deg)

Kθ (kJ mol−1)

O−S−C S−C−H H−C−H

110 110 110

503 450 503



RESULTS AND DISCUSSION Thermodynamic and Structural Properties. The bond lengths and the bending and torsional angles computed from our ab initio quantum calculation, which are gathered in Tables 1, 2, and 3, are in good agreement with the structural parameters obtained by the groups of Durig,7 Wang,8 and Carvalho,9 using similar ab initio procedures at several levels and with different basis sets. According to the work of Durig et al.,7 we have considered the gauche conformation of the OH group, which corresponds to a local minimum, in which the hydrogen atom is close to eclipsing one of the oxygen atoms. A plot of the most favorable structure is depicted in Figure 1. As it has also been reported by these authors,7−9 the conformation of the CH3 group is staggered. The MD density value at 298 K and 1 atm is 1.462 ± 0.004 g/cm3, which is in a very good agreement with the experimental data34 at the same temperature and pressure 1.475 g/cm3. The intermolecular structure of liquids is normally described in terms of the center of mass-center of mass gCM−CM(r) and the partial radial distribution functions. In the case of MSA the most interesting correspond to the hydrogen−oxygen and oxygen−oxygen pairs (gH−O(r) and gO−O(r), respectively). The gCM−CM(r) functions calculated at 200, 298, and 400 K have been plotted in Figure 2 (top). Three clear peaks, located at 0.50, 0.95, and 1.35 nm, can be observed. The temperature

Table 3. Dihedral Potential Parameters (Multiplicity 3) dihedral

angle (deg)

Kφ (kJ mol−1)

H−C−S−O,O′ H′−O′−S−O,C

0 145

2.93 450

nonbonded van der Waals interactions have been calculated using the Lennard-Jones 12−6 potential. The OPLS-AA homoatomic pair parameters,26 which are collected in Table 4, have been employed. The heteroatomic pair parameters are calculated from the homoatomic ones using the well-known Lorentz−Berthelot mixing rules. The electrostatic interactions Table 4. Lennard-Jones Parameters and Partial Charges atom

σ (nm)

ε (kJ mol−1)

q (e)

H C S O O′ H′

0.250 0.350 0.355 0.296 0.312 0

0.1255 0.2761 1.0460 0.8786 0.7113 0

0.06 −0.18 0.84 −0.40 −0.48 0.44 3424

dx.doi.org/10.1021/jp500817s | J. Phys. Chem. B 2014, 118, 3423−3430

The Journal of Physical Chemistry B

Article

Figure 1. Most favorable structure of methanesulfonic acid. Labels correspond to the atom types defined in Table 4.

Figure 3. Hydrogen−oxygen radial distribution function gH−O(r) at 298 K and hydroxyl and sulfonyl contributions (top). gH−O(r) at 200, 298, and 400 K (bottom).

bonding distance r(H−O) = 0.1784 nm predicted by Durig et al.7 for a MSA dimer using ab initio calculations. Moreover, at short distances (

Molecular dynamics simulation study of methanesulfonic acid.

A molecular dynamics simulation study of methanesulfonic acid has been carried out using a reliable force field in a large range of temperatures. Seve...
1MB Sizes 0 Downloads 3 Views