Molecular Conformations and the Liquid Structure in Bis(methylthio

Jun 29, 2012 - The structure of PX 3 (X = Cl, Br, I) molecular liquids from X-ray diffraction, molecular dynamics simulations, and reverse Monte Carlo...
0 downloads 0 Views 2MB Size
Download PDF
Article pubs.acs.org/JPCB

Molecular Conformations and the Liquid Structure in Bis(methylthio)methane and Diethyl Sulfide: Diffraction Experiments vs Molecular Dynamics Simulations Orsolya Gereben and László Pusztai* Institute for Solid State Physics and Optics, Wigner Research Centre for Physics, Hungarian Academy of Sciences, P.O.B 49, H-1525 Budapest, Hungary ABSTRACT: Molecular dynamics computer simulations with the OPLS-AA and EncadS all-atom force fields have been performed for liquid bis(methylthio)methane and diethyl sulfide. For bis(methylthio)methane, additional simulations using gas electron diffraction and ab initio bond length and bond angle parameters were also performed. Although there were small differences between the structures provided by the different molecular dynamics simulations of bis(methylthio)methane, the main conformer was found to be the AG one (54%), followed by the G+G+ one (31%) in all the calculations. This is in contrast with gas phase electron diffraction and ab initio calculations, where the main conformer was G+G+ (70%; beside that, 30% AG was assumed). For diethyl sulfide, no experimental data was present for the gas phase and the different interatomic potential sets produced different conformer ratios: applying OPLS-AA resulted in 51% AG and 40% AA, while from EncadS the main conformer was G+G+, with an occurrence of 69.5% (AG was present with 29.6%). Comparing the MD calculated structure factors to results of X-ray diffraction experiments, some differences could be found at Q < 5 Å−1 in all cases. For bis(methylthio)methane, the EncadS potential produced the X-ray structure factor with the smallest difference to the experimental one, whereas for diethyl sulfide the OPLS-AA parameter set proved to be the more successful. The EncadS force field provided a better approximation to the experimental value of the enthalpy of vaporization for bis(methylthio)methane; for diethyl sulfide, the two different force fields gave similar results.

1. INTRODUCTION Diethyl sulfide [DEtS, IUPAC name: (ethylthio)ethane] and bis(methylthio)methane [BMtMe, IUPAC name: 2,4-dithiapentane] are both compounds with distinctive odor. BMtMe’s odor resembles the smell of fungi, and is applied as food flavor for mushroom, tomato, and truffle.1 DEtS can be found in cabbage, broccoli, and cauliflower, and it has a garlic smell; this material is also used in plating baths for gold and silver and as a solvent for anhydrous minerals.2 The gas phase structure of BMtMe was investigated by Page3 by gas phase electron diffraction (GED) and ab initio calculations. According to their findings, BMtMe in the gas phase can be best described as a mixture of 70% G+G+ and 30% AG conformers (for the definition of conformers, see Figure 1, ref 3). For obtaining experimental information on the liquid structure, high-energy X-ray diffraction (HEXRD) measurements have been conducted. Due to the high hydrogen content of these materials, neutron diffraction is not yet practical for gathering extra amounts of direct information on the structure. This work reports on the results of molecular dynamics (MD) simulations for both materials in their liquid phase; comparisons of these simulation results with each other as well as with the experimental structure factor are also provided. The original goals were (a) the calculation of the structure factor from the MD trajectories and (b) generating suitable starting © 2012 American Chemical Society

Figure 1. Schematic drawings of the anti and gauche conformers, shown on the anti−anti conformer as an example. The backbone of the molecule is represented by circles; only the relevant ligands are displayed.

configurations for subsequent structural modeling studies; during the course of the investigations, however, the spotlight has been shifted to the conformational analyses of the molecules in the two liquids. The present contribution therefore focuses on the agreement between simulation and experiment, as well as on the different molecular conformations, as a function of the actual force field applied. In the near future, reverse Monte Carlo (RMC) refinement4−6 of the Received: April 25, 2012 Revised: June 6, 2012 Published: June 29, 2012 9114

dx.doi.org/10.1021/jp3040117 | J. Phys. Chem. B 2012, 116, 9114−9121

The Journal of Physical Chemistry B

Article

Table 1. Lennard-Jones Parameters and Partial Charges for the Atom Types Used in the MD Simulations for the Two Liquids OPLS-AA

EncadS

name

σ (Å)

ε (kJ/mol)

Q (e)

σ (Å)

ε (kJ/mol)

q (e)

C (CH3,) C (CH2) S (BMtMe) H

3.5 3.5 3.6 2.5

0.276144 0.276144 1.48532 0.12552

−0.013 (BMtMe), −0.18 (DEtS) 0.216 (BMtMe), 0.048 (DEtS) −0.335 (BMtMe), −0.336 (DEtS) 0.06

3.84423 3.84423 3.84423 2.54129

0.308863 0.308863 0.308863 0.158992

−0.357 −0.238 0.0 0.119

structure will be conducted, so that structural models that are fully consistent with X-ray diffraction data may be analyzed.

Δhvap(T , p) = ugas(T ) − uliq(T , p) + RT

(1)

where Δhvap is the enthalpy of vaporization, ⟨Δugas⟩ and ⟨Δuliq⟩ are the molar internal energies of the ideal gas and the neat liquid, and R is the gas constant. The comparison of the experimental and MD-based F(Q) was performed according to

2. SIMULATION DETAILS Molecular dynamics simulations have been performed by the GROMACS simulation package (version 4.0)7 at T = 293 K in the NVT ensemble. Two different all-atom force fields, named OPLS-AA8 and EncadS,9 operating with full solvent charges, were used for both liquids. In each calculation, 1331 molecules were put in cubic simulation cells, with atomic number densities of 8.38345022 × 10−2 and 7.6605545 × 10−2 Å−3 which corresponds to (experimentally determined) bulk densities of 0.837 and 1.059 g/cm3 for DEtS and BMtMe, respectively. Preliminary simulations were performed by keeping the temperature constant by either the Berendsen thermostat10 with temperature coupling time constant τT = 0.1 ps or by the Nosé− Hoover thermostat11 with τT = 0.5 ps. Time steps of 1 and 2 fs were both tested. As there were no differences regarding the outcome of the simulation, structures resulting from the 2 fs time step simulation, using the Berendsen thermostat, will be presented. Both the OPLS-AA8 and EncadS9 force fields use pair interactions that are composed of (a) 6−12 Lennard-Jones (LJ) interactions for representing dispersion and repulsion effects and (b) Coulomb interactions. LJ parameters and partial charges are given in Table 1. Cutoffs for the LJ and for the short-range Coulomb interactions were set at 11 and 9 Å, respectively. For the long-range part of the Coulomb interaction, the particle-mesh Ewald12 summation, with interpolation order 4 and a grid size of 0.12 nm, was applied. The relative strength of the electrostatic interaction at the cutoff was ∼10−5, providing an accuracy of 5 × 10−3 which is still higher than the accuracy of the Lennard-Jones calculations. The leapfrog algorithm was used for integrating Newton’s equations of motion with a time step of dt = 2 fs, as mentioned above. In each case, the total energy reached its equilibrium value within 100 ps. MD configurations were collected 20 ps apart, and 76 time frames from the last 1500 ps of the 2000 ps trajectory were used for calculating the X-ray weighted structure factor, F(Q), and the bond length, bond angle, and dihedral angle distributions, if not stated otherwise. The histogram calculation and the Fourier transformation was performed exactly the same way as in our reverse Monte Carlo code,13 after modifying the g_rdf software of the GROMACS7 package for this purpose. The enthalpy of evaporation was calculated at the boiling point of both materials, which is 420 K for BMtMe and 365 K for DEtS, at 105 Pa (ref 14). First, a liquid phase simulation for the whole system was performed for 200 ps and then a gas phase simulation for one molecule for 20 ps without periodic boundary conditions. The enthalpy of vaporization was calculated on the basis of the work of Kelkar15 according to the formula

χ2 =

Npi

∑ (aFjE + b + cQ j + dQ j 2 − F Cj )2 j=1

(2)

where a, b, c, and d are renormalization constants. If renormalization was applied to the experimental structure factor, then parameters a, b, c, and d were determined to give the minimum of χ2; if no renormalization was applied, then a = 1, b = 0, c = 0, and d = 0 were used. All MD simulations used flexible molecules (i.e., none of the bond lengths or bond angles were fixed): harmonic bond length and bond angle terms were used. Dihedral angles were calculated as Ryckaert−Bellemans dihedrals,16 according to the GROMACS implementation of the OPLS8 force field. For EncadS,9 proper dihedrals were applied, as described in the GROMACS manual.7 For the determination of the equilibrium liquid structure of BMtMe, four MD simulations were carried out. The first and second calculations applied the original parameters of the two force fields; these will be denoted as simulations “OPLS-AA” and “EncadS”. For the third and fourth simulations, parameters were taken from Page’s gas phase electron diffraction and ab initio study;3 these calculations will be referred to as “OPLSAA_Page” and “EncadS_Page”. For these latter cases, force constants and any other missing parameters were taken from the original OPLS-AA and EncadS force fields, respectively. No SCS angular potential parameters were provided in either the OPLS-AA or the EncadS force field; as a first approximation, force constants for the CCS angle of the given force field were used, along with the experimental value of the bond angle. For DEtS, no applicable experimental data could be found; therefore, only two MD simulations, “OPLS-AA” and “EncadS”, were performed. The bond length and bond angle parameters are listed in Table 2. Dihedral angles were used as they are defined in the OPLSAA and EncadS force field files, even for simulations that otherwise take the experimental bond length and bond angle values. As no value exists in the OPLS-AA force field for the CSCS dihedral in BMtMe, the CCSC dihedral angle value was taken instead. The bin size for the cosine distributions of bond angles was d(cos θ) = 0.05; in the dihedral angle distributions, bins were dψ = 1° wide. The value halfway between the average gauche and average anti peak position of the dihedral distributions served as a limiting value when classifying a given conformer as gauche or anti. 9115

dx.doi.org/10.1021/jp3040117 | J. Phys. Chem. B 2012, 116, 9114−9121

The Journal of Physical Chemistry B

Article

Table 2. Equilibrium Bond Length and Bond Angle Parameters for the Two Liquids as Applied in the MD Calculations bond CH CH CH CS CS CS (methyl) CS (methylene) CC CC angle HCH HCH HCH HCS HCS HCS CSC CSC CSC SCS HCC HCC CCS CCS

be (Å) 1.09 1.09 1.108 1.81 1.808 1.805 1.806 1.529 1.525 θ (deg) 107.8 109.5 108.9 109.5 110.0 107.5 98.9 99.1 102.8 115.9 110.7 109.5 114.7 115.5

source

Figure 2. Labeling of atoms in the BMtMe and DEtS molecules.

OPLS-AA EncadS Page OPLS-AA EncadS Page Page OPLS-AA EncadS source

partial charge distributions of the different force fields (Table 1). The enthalpy of vaporization was calculated by using the original OPLS-AA and EncadS force fields, giving Δhvap(420 K,105 Pa) = 52.3 ± 0.3 kJ/mol for the OPLS-AA and 35.4 ± 0.4 kJ/mol for the EncadS force field. The experimental value is 36.816 kJ/mol,20 so the value of the EncadS force field was able to reproduce the experimental value much better. Next we compare the various equilibrium average molecular structures produced by the MD simulations to each other, as well as to the GED/ab initio structure of Page.3 It has to be kept in mind that the molecular structure may be different in the gaseous and liquid phases. Considering the equilibrium bond lengths of the different MD simulations, there are no large variations. The two MD using the original force field parameters are closer to each other than the two “MD_Page” simulations: the former display slightly smaller CH (1.09 Å) and larger CS (1.81 Å) bond distance values. On the other hand, the “...-Page” calculations produce a better agreement with experimental data (see Table 3). The average MD HCH and SCS angles are in good agreement with the experimental values: the HCS MD angles are all a bit larger than the corresponding experimental ones. The largest difference is found for the CSC angle (Figure 4a), where the “OPLS-AA” average angle is 3.7° smaller; this value will increase to almost the experimental value in the “OPLS-AA_Page” simulation. The EncadS potential brings about an ∼1.5° smaller angle, whereas the “EncadS_Page” run, a 1.6° larger CSC angle; the shape of the distributions are very similar. For all the other angle distributions, the EncadS force field produces sharper distributions than the OPLS-AA one (not shown). The HCSC MD dihedral distribution and average dihedral angles for the gauche and anti forms are very similar, close to 60 and 180°; unfortunately, no experimental or ab initio data are available for comparison. For the CSCS dihedral angle, the OPLS-AA force field produced a larger average gauche dihedral angle (65−66°) than the EncadS potential (60°). The 54° of Page,3 approximated from a least-squares fit to the GED data, is much lower. The average dihedral angle of ∼69° would be expected from the ab initio calculation on the 70% G+G+-30% AG mixture.3 The A/G conformer ratio does not differ much for the two MD simulations (40/60% for “OPLS-AA” and 35/ 65% for “EncadS”), although distributions from the “EncadS” run are sharper (Figure 4b). On the basis of the GED-ab initio data, this ratio comes out to 15/85%. In all the MD studies, the main conformer is AG (∼54%), the second is G+G+(G−G−) (31−39%), followed by AA (8−13%) and a very small amount (if any at all) of G+G− (0−4%). The GED-ab initio study suggested that G+G+(G−G−) is the main conformer with an occurrence of 70%, whereas 30% AG is supposed to be present (to achieve the best fit to measured GED data). The stability of the G+G+(G−G−) form is probably the result of the anomer effect, which is caused by the delocalization of a lone electron pair of the heteroatom X to the σ*(C−X) orbital of the neighboring C atom.21 Another explanation may be that in the anti conformer the dipoles created by the polarization of

OPLS-AA EncadS Page OPLS-AA EncadS Page OPLS-AA EncadS Page Page OPLS-AA EncadS OPLS-AA EncadS

Detailed conformer analysis was performed for the last configurations of the simulations, but the calculation was performed for another, 400 ps earlier configuration of the equilibrium, as well. The difference between the conformer ratios of the two configurations was less than 1%; i.e., the systems have certainly been well equilibrated. Concerning experimental data, high energy X-ray diffraction measurements were carried out at the SPring-8 synchrotron radiation facility (Hyogo, Japan), using the single-detector diffractometer setup of the BL04B2 (high-energy X-ray diffraction) beamline,17 as parts of a series of measurements on sulfur-containing molecular liquids.18 Corrections to yield structure factors have been made by standard procedures.19 More details can be found in ref 18.

3. RESULTS AND DISCUSSION 3.1. Bis(methylthio)methane, BMtMe. For BMtMe, four MD simulations have been performed (see section 2). For the two pairs of simulations (OPLS-AA−OPLS-AA_Page; EncadS−EncadS_Page), the average energies were almost identical within the pairs (21.9, 23.1, 125.8, and 124.4 kJ/ mol, in the above order). Note that “OPLS-AA” simulations produced much lower energies; the majority of the difference comes from the reciprocal space Coulomb contribution. Detailed analysis was performed to see whether the difference could be attributed to any specific groups of atoms. It could be identified that the large increase of the reciprocal space Coulomb energy in simulation “EncadS” was caused by the CAH3−CBH2, S−S, and intermolecular CBH2 interactions (see labeling in Figure 2); these were somewhat moderated by the decrease of the CAH3−S and CBH2−S contributions. That is, the large energy difference can be attributed to the different 9116

dx.doi.org/10.1021/jp3040117 | J. Phys. Chem. B 2012, 116, 9114−9121

The Journal of Physical Chemistry B

Article

Table 3. BMtMe Equilibrium Structural Parameters for the Different MD Simulations and as Obtained from Experimental Electron Diffraction and ab Initio Calculation Data of Page OPLS-AA

OPLS-AA_Page

EncadS

CH (Å)

1.090

1.108

1.091

1.109

CS (Å) HCH (deg) HCS (deg) CSC (deg) SCS (deg) HCSC (deg) CSCS (deg)

1.810 108.1 109.7 99.1 114.5 ±59, 180 ±66 (60.5%), 180 (39.5%) 31.3 53.9 12.8 1.9

1.804 109.9 108.4 102.7 115.3 ±58, 180 ±65 (61.7%), 180 (38.4%) 31.4 53.1 11.8 3.7

1.810 108.8 109.3 101.3 114.3 ±60, 180 ±60 (64.8%), 180 (35.2%) 37.2 55.2 7.6 0.0

1.807 109.8 108.4 104.4 116.1 ±60.5, 180 ±60 (65.7%), 180 (34.3%) 38.8 53.0 7.7 0.4

G+G+, G−G− (%) AG AA G+G−

Page3

EncadS_Page

1.1080(CH3), 1.1098(CH2) 1.805(CH3S), 1.806(CH2S) 108.9 107.5 102.8 115.9 ±54 (85%), 180 (15%) 70 30

Figure 3. The AG, G+G+, AA, and G+G− conformers of BMtMe. The first three pictures are from the final configurations of calculation “EncadS”; the last picture is from simulation “EncadS_Page”.

Figure 4. (a) CSC angle and (b) CSCC dihedral angle distributions as obtained from the various MD simulations of BMtMe. In part a, the peaks are following each other from left to right in the order of the legends. In part b, the two OPLS-AA distributions run together (smaller peaks), and the higher peaks are the two EncadS distributions that run together.

the X and Y atoms of the C−X−C−Y part are aligned at the same direction, therefore repelling each other, while in the gauche form they are turned almost opposite to each other and this way they are capable of stabilizing this, stereochemically otherwise unfavorable, position. It is not reasonable to expect that MD simulations without polarizable potentials could reproduce this effect. Marked differences can be spotted in terms of the nonbonding CA−CA_nb distance distributions (Figure 5), not only between results produced by different force fields but also between the original and “_Page” simulations of the same force field. Here, for the AG and AA conformers, the peaks are shifted toward larger r values for the “_Page” simulations. This latter effect can be attributed to the larger SCS and CSC angles of the “_Page” simulations, as compared to their counterparts that use the original force field parameters. The “EncadS”

Figure 5. The nonbonding CA−CA_nb distance distribution in BMtMe, as obtained from applying different potential parameters.

simulation produced well-defined peaks that can be explained by the sharper dihedral distributions brought about by this force field. CSCS dihedral angles were calculated separately for the different conformers found in the final configuration of each 9117

dx.doi.org/10.1021/jp3040117 | J. Phys. Chem. B 2012, 116, 9114−9121

The Journal of Physical Chemistry B

Article

Table 4. Average Dihedral Angles for the Different Conformers Found in the Final Configurations of the Various MD Simulations of Bis(methylthio)methane (BMtMe) and Diethyl Sulfide (DEtS)a AA

A (AG)

OPLS-AA (66) OPLS-AA_Page (65) EncadS (60) EncadS_Page (60) HF/6-311 + G(d)3 MP2/6-311 + G(d)3

181.2(18.7) 177.6(16.9) 180.5(8.8) 180.0(9.4) 179.54 158.59

179.8(18.1) 180.9(18.8) 180.1(8.7) 179.5(9.0) 186.48 191.48

OPLS-AA (77) EncadS (60.5)

179.9(18.4) 183.6(7.7)

179.2(19.5) 180.9(9.2)

G (AG) BMtMe, CSCS 69.2(15.0) 67.2(16.1) 60.4(8.3) 59.2(8.8) 74.14 74.8 DEtS, CCSC 78.8(15.3) 62.8(9.5)

G+G+, G−G−

G+G−

68.5(14.1) 67.5(16.1) 60.8(8.3) 59.4(8.4) 67.52 67.69

−94.6(20.6), 82.6(22.6), 69.6(13.0), 107.6(10.7) −88.5(21.3), 85.6(20.7), 69.5(11.6), 104.5(11.9) −70.5(12.3), 73.0(9.1), 66.6(9.9), 76.8(9.3) −83.83, 83.83 −82.84, 82.84

79.9(14.8) 60.9(7.2)

−85.3(16.8), 86.0(16.8), 73.1(12.4), 98.2(9.8) −96.4(9.3), 60.5(4.6), 60.5(4.6), 96.4(9.3)

a

Angles are shown in degrees; values in round brackets are the corresponding standard deviations. Average gauche angles are given in round brackets following the simulation name. For the G+G− conformation, first the separate minus and plus dihedral angle averages, then the smaller and larger absolute value dihedral angle averages (in italics) are displayed.

Figure 6. X-ray weighted structure factors of BMtMe. (a) F(Q)'s for the four MD simulations as compared to the experimental data. (b) The original and renormalized experimental F(Q) as compared to the same function from simulation “EncadS”.

according to eq 2 (with and without renormalization) are given in Table 5. Visually, there is not much difference between the

MD calculation; the purpose was to see whether the conformers have any preference for a given section of dihedral distribution curve (Table 4). The average gauche angle for the simulation is also given; the average dihedral angle was 180° in every case (see Tables 3 and 6). In the case of the conformers where only one gauche angle can be found in the molecule [AG, G+G+(G−G−)], the absolute values of the angles were averaged, as there was hardly any difference between the separate averages for the “minus” and “plus” angles. For the G+G− conformer, the minus and plus angles were averaged separately. As there was sometimes a substantial difference between the minus and plus gauche angle averages, the smaller and larger absolute value angles for every molecule were averaged separately as well. For the anti conformation, the average dihedral angle for both the AA and AG conformers was close to the average 180° value. Ab initio calculations3 for the gas phase suggest a more distorted structure. The dihedral angle in the gauche conformation AG and the two such angles in G+G+(G−G−) are close to each other for a given simulation, as well as to the average gauche angle of the simulation, particularly for the EncadS force field. Analogue values for the OPLS-AA force field are closer to the ab initio values. The G+G− conformer is an energetically less favorable one, according to both the ab initio3 and our MD calculations; it is also the farthest from the average gauche angle, indicating a distorted gauche−gauche structure. The MD calculation shows that there is a preference for two different, one smaller and one larger, gauche dihedrals within a molecule, while the ab initio calculation suggested two equal angles. The average X-ray total scattering structure factors calculated for all four MD simulations (without any renormalization) are shown in Figure 6a; the corresponding χ2 values calculated

Table 5. Squared Differences, χ2, between Average F(Q)'s Resulting from Several MD Simulations of BMtMe and the Experimental Data (without and with Renormalization) simulations

χ2

χ2 renorm

OPLS-AA OPLS-AA_Page EncadS EncadS_Page

20.0 22.6 18.8 21.3

13.0 15.1 9.1 11.8

MD structure factors; their differences as compared to the experimental results are much larger. The EncadS F(Q) is closest to the experimental data, especially when using the renormalization option (χ2 = 9.1). Figure 6b shows F(Q) from the EncadS simulation (with and without renormalization) in comparison with the experimental data. The renormalization decreased the difference between simulation and experiment in the regions before and around the first maximum but not substantially so around the second maximum and between the first and second peaks. Simulations using the GED-ab initio data produced worse agreement with the experimental F(Q) than simulations using the original force field for both materials. The same tendency was observed in our earlier study of dimethyl sulfide, dimethyl disulfide, and dimethyl trisulfide.18 Comparing partial radial distribution functions (prdf's) of the four simulations (Figure 7), it may be concluded that beyond the first maximum both the “OPLS-AA” and “OPLS-AA_Page” and the “EncadS” and “EncadS_Page” curve pairs run together. This is in agreement with the expectations that the larger distance region is not influenced by the exact (and not very 9118

dx.doi.org/10.1021/jp3040117 | J. Phys. Chem. B 2012, 116, 9114−9121

The Journal of Physical Chemistry B

Article

Table 6. Equilibrium Bond Length, Bond Angle, and Dihedral Angle Values for the Different MD Simulations of DEtSa OPLS-AA CC (Å) CH (Å) CS (Å) HCH (deg) HCC (deg) HCS (deg) CSC (deg) CCS (deg) HCSC (deg) CSCS (deg)

1.53 1.09 1.81 107.51 110.23 108.52 100.07 113.98 ±60.5, 180 ±77 (34.4%), 180 (65.6%)

G+G+, G−G− (%) AG (%) AA (%) G+G− (%)

5.9 51.0 40.0 3.1

a

EncadS 1.52 1.09 1.81 109.14 108.84 108.79 102.44 114.59 ±60.5, 180 ±60.5 (84.5%), 180 (15.5%) 69.5 29.6 0.7 0.2

The occurrences of the various conformers are also shown.

Figure 7. Partial radial distribution functions of BMtMe.

different) bond length and bond angle values. On the other hand, around the first maximum of the C−C, C−H, and H−H prdf's, the two “MD_Page” simulations are closer to each other than either the two “OPLS” or the two “EncadS” ones. The CS and CH prdf's are virtually identical for the four simulations, and there are no significant differences in terms of the S−H and H−H prdf's, either. There is a pronounced second peak of the C−C prdf on the “EncadS” simulation result that may be attributed to the presence of the G+G+ conformer (for details, see Figure 5). The largest differences can be found in the intermolecular part of the S−S prdf: the “EncadS” simulations produce a peak at ∼3.85 Å, indicating a tighter arrangement of molecules. 3.2. Diethyl Sulfide, DEtS. For DEtS, two MD simulations, one using the OPLS-AA while the other using the EncadS force field, have been carried out. The EncadS force field produced a structure characterized by a deeper energy value: the average total energy is 78.6 and 26.5 kJ/mol for the OPLS-AA and EncadS force fields, respectively. Detailed analysis revealed that the difference originates from the reciprocal space Coulomb contribution, and mainly from the 1−4 CH3, CH2 groups. Thus, the difference in terms of the energy is caused by the different charge distributions, similarly to the case of BMtMe (see Table 1 for the partial charges). The enthalpy of vaporization was calculated at the boiling point for both force fields, providing Δhvap(365 K, 105 Pa) = 36.0 ± 0.3 kJ/mol for OPLS-AA and 35.4 ± 0.3 kJ/mol for the EncadS force field. The experimental value of 31.77 kJ/mol was reported in ref 22. Both values are close to the experimental value, the value of the EncadS force field being somewhat more favorable. In terms of the equilibrium bond length values, the two simulations produced similar results (see Table 6). Larger differences can be found for the equilibrium angles, the CSC (Figure 8a) and HCH angle differences being the largest by 2.37 and 1.63°, in favor of “EncadS”. Remarkable differences can be observed for the CSCS dihedrals: the gauche “OPLSAA” dihedral angle is much larger, ±77°, than its “EncadS” counterpart, ±60.5° (Figure 8b). It is also visible that the ratios

Figure 8. CSC angle and CCSC dihedral angle distributions for DEtS, as obtained from MD simulations using the OPLS-AA and EncadS force fields.

of the gauche and anti 1−4 C−C atoms are very different: the gauche conformation (84.5%) is preferred by “EncadS”, whereas the anti one (65.6%) by “OPLS-AA”. The most interesting difference between the two equilibrium structures is the remarkably different ratios of the four possible conformers. In this case, contrary to BMtMe, even the order of the most preferable conformers is different: while in the case of the OPLS-AA force field there are two main conformers, the AG and the AA with 51 and 40%, respectively, and small amounts, 5.9 and 3.1% of G+G+(G−G−) and G+G−, in the case of the EncadS force field, the dominant conformer is G+G+(G−G−) with 69.5%, while 29.6% of AG and hardly any (0.7 and 0.3%) AA and G+G− can be found (Table 6). An example of the conformers can be seen in Figure 9. The average CCSC dihedral angles in the final configuration of the two MD simulations were calculated separately for the conformers, similarly to liquid BMtMe (Table 4). The pattern is very similar to that found for BMtMe; namely, the anti conformers in AA and AG are close to the average anti angle of 180°, whereas the gauche angle(s) in AG and G+G+(G−G−) are close to the average gauche dihedral angle of the given simulation. The stereochemically, and therefore energetically, less favorable G+G− conformer contains highly distorted gauche angles. The X-ray structure factors calculated from the MD configurations reproduce the main features of the experimental data quite well (see Figure 10a), although there are visible differences especially around the first, second, and third maxima. As it is shown in Table 7, the F(Q) calculated from 9119

dx.doi.org/10.1021/jp3040117 | J. Phys. Chem. B 2012, 116, 9114−9121

The Journal of Physical Chemistry B

Article

Figure 9. The AG, AA, G+G+, and G+G− conformers as found in the final configurations of the MD simulation using the OPLS-AA potential set for DEtS.

Figure 10. X-ray weighted structure factors of DEtS. (a) F(Q)'s calculated from the MD simulations as compared to the experimental results. (b) F(Q) calculated from the OPLS-AA simulation as compared with the experimental and renormalized experimental F(Q).

Table 7. Squared Differences, χ2, between F(Q)'s Resulting from the Two MD Simulations of DEtS and the Experimental Data (without and with Renormalization) simulations

χ2

χ2 renorm

OPLS-AA EncadS

10.5 20.1

5.9 10.5

the “OPLS-AA” simulation is closer to the experimental data than that derived from run “EncadS”. If renormalization was applied, according to eq 2, then χ2 decreased to about half of its original value. Figure 10b shows the best fit resulting from the “OPLS-AA” simulation, the experimental and renormalized experimental data. As it is visible, renormalization mainly decreased the difference before and around the first maximum but had no effect on the third peak (around 5 Å−1). Comparing prdf's of the two MD simulations, we can conclude from Figure 11 that differences occur either in the intermolecular or in the larger intramolecular distance regions of the molecule but not in the region of bonding distances. The C−S and C−H partials are almost identical, similar to the case of BMtMe. There is an extra peak in the C−C prdf at 3.3 Å for simulation “EncadS”, which results from the gauche CA−CB interaction (see labeling in Figure 2). On the same prdf, around 4.2 Å, simulation “OPLS-AA” produced a slightly larger peak, coming solely from the anti CA−CB interaction; the lower maximum of the curve corresponding to calculation “EncadS” in the same region can be attributed to the anti CA−CB and G+G+ CA−CA interactions, as can be seen in Figure 12. The intermolecular S−S peak from simulation “OPLS-AA” is much more pronounced, indicating a more uniform alignment of the molecules. The maximum of the H−H partial at 4 Å for calculation “EncadS” results from the hydrogen atoms of the gauche CA−CB arrangement, as detailed analysis of the H−H nonbonding distance distribution revealed (Figure 12).

Figure 11. Partial radial distribution functions of DEtS, as calculated from results of MD simulations using the OPLS-AA and EncadS force fields.

Figure 12. (a) CA−CB and (b) CA−CA nonbonding distance distributions in liquid DEtS, as calculated from MD simulation results using the OPLS-AA and EncadS force fields.

4. SUMMARY AND CONCLUSIONS All-atom molecular dynamics studies have been performed for bis(methyltio)methane (BMtMe) and diethyl sulfide (DEtS) in 9120

dx.doi.org/10.1021/jp3040117 | J. Phys. Chem. B 2012, 116, 9114−9121

The Journal of Physical Chemistry B

Article

(11) (a) Nose, S. Mol. Phys. 1984, 52, 255. (b) Hoover, W. G. Phys. Rev. A 1985, 31, 1695−1697. (12) Essman, U.; Perela, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, K. G. J. Chem. Phys. 1995, 103, 8577−8592. (13) Gereben, O.; Jóvári, P.; Temleitner, L.; Pusztai, L. J. Optoelectron. Adv. Mater. 2007, 9, 3021−3027. (14) National Institute of standard and technology, http://webbook. nist.gov/chemistry/. (15) Kelkar, M. S.; Manning, E. J. J. Phys. Chem. B 2007, 111, 9424− 9427. (16) Ryckaert, J. P.; Bellemans, A. Faraday Discuss. Chem. Soc. 1978, 66, 95−99. (17) Kohara, S.; Itou, M.; Suzuya, K.; Inamura, Y.; Sakurai, Y.; Ohishi, Y.; Takata, M. J. Phys.: Condens. Matter 2007, 19, 506101− 506115. (18) Gereben, O.; Kohara, S.; Pusztai, L J. Mol. Liq. 2012, 169, 63− 73. (19) Kohara, S.; Suzuya, K.; Kashihara, Y.; Matsumoto, N.; Umesaki, N.; Sakai, I. Nucl. Instrum. Methods Phys. Res., Sect. A 2001, 467−468, 1030−1033. (20) LookChem, http://www.lookchem.com/Bis-methylthiomethane/. (21) (a) Kirby, A. J. The Anomeric Effect and Related Stereoelectronic Effects at Oxygen; Springer-Verlag: Berlin, 1983. (b) Kneisler, J. R.; Allinger, N. L. J. Comput. Chem. 1996, 17, 757−766. (c) Smith, G. D.; Jaffe, R. L.; Yoon, D. Y. J. Phys. Chem. 1994, 98, 9072−9077. (22) Enthalpies of Vaporization of Organic Compounds: A Critical Review and Data Compilation; Blockwell Scientific Publications: Oxford, U.K., 1985.

the liquid phase, applying the OPLS-AA and EncadS force fields. For BMtMe, two other simulations using the GED and ab initio bond length and bond angle data of Page3 were also performed. For BMtMe, the OPLS-AA, while for DEtS, the EncadS potential parameter sets produced lower energy structures. Comparing X-ray weighted F(Q) of the simulations to the corresponding experimental data, for BMtMe calculation “EncadS” using the original force field parameters, whereas for DEtS simulation “OPLS-AA” yielded the best agreement. The enthalpy of vaporization was also calculated; in the case of BMtMe, the value obtained for the EncadS force field was much closer to the experimental one, similarly to the X-ray structure factor results. In the case of DEtS, the two calculated values were close to each other, so in this respect no substantial difference can be found between the force fields. Detailed conformation analyses revealed that, although there were differences between the structures produced by the different force fields for BMtMe, the main conformer was AG (54%), followed by G+G+ (31%) in all the simulations; in the gas phase, on the other hand, the main conformer is G+G+ (70%) and 30% AG was assumed. As there had been no gas phase experimental data available for DEtS, only two simulations were performed which produced markedly different conformer ratios: calculation “OPLS-AA” resulted in 51% AG and 40% AA, while in calculation “EncadS” the main conformer was G+G+ (69.5%) followed by AG (29.6%). On the basis of the above findings, none of the two force fields can be considered to be superior for describing the structure of the two liquids in question; for a definite statement, better agreement with experimental data is required. This is the main reason why detailed reverse Monte Carlo modeling investigations are underway.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We thank Dr. S. Kohara (Spring-8, Japan) for his invaluable help with the X-ray diffraction experiments. Financial support was provided by the Hungarian National Basic Research Fund (OTKA), via Grant No. 83529.

■

REFERENCES

(1) http://www.thegoodscentscompany.com/data/rw1057521. html#. (2) Merck index of Chemicals and Drugs, 9th ed., monograph 3801. (3) Page, E. M.; Rice, D. A.; Aarset, K.; Hagen, K.; Genge, A. R. J. J. Phys. Chem. A 2000, 104, 6672−6676. (4) McGreevy, R. L.; Pusztai, L. Mol. Simul. 1988, 1, 359−367. (5) McGreevy, R. L. J. Phys.: Condens. Matter 2001, 13, R877−R914. (6) Evrard, G.; Pusztai, L. J. Phys: Condens. Matter 2005, 17, S1−S13. (7) van der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. C. J. Comput. Chem. 2005, 26, 1701−1718 ( http://www.gromacs.org). (8) Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. J. Am. Chem. Soc. 1996, 118, 11225−11236. (9) (a) Levitt, M.; Hirshberg, M.; Sharon, R.; Daggett, V. Comput. Phys. Commun. 1995, 91, 215−231. (b) Levitt, M. J. Mol. Biol. 1983, 168, 595−620. (10) Berendsen, H. J. C.; Postma, J. P. M.; DiNola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684−3690. 9121

dx.doi.org/10.1021/jp3040117 | J. Phys. Chem. B 2012, 116, 9114−9121