Intermolecular potential functions and Monte Carlo simulations for

Purdue University, West Lafayette, Indiana 47907 (Received: May 27, 1986;. In Final Form: July 24, ... new parameters which were chosen to yield g...
2 downloads 0 Views 1MB Size
J. Phys. Chem. 1986, 90, 6379-6388

6379

Intermolecular Potential Functlons and Monte Carlo Simulations for Liquid Sulfur Compounds William L. Jorgensen Department of Chemistry, Purdue University, West Lafayette, Indiana 47907 (Received: May 27, 1986; In Final Form: July 24, 1986)

Intermolecular potential functions have been developed for use in computer simulations of liquid thiols, sulfides, and disulfides as well as the sulfur-containing side chains of proteins. The simple OPLS model was used and required introduction of few new parameters which were chosen to yield good descriptions of the structures and energetics for bimolecular and ion-molecule complexes and for the thermodynamic properties of liquid sulfur compounds. The principal testing involved Monte Carlo statistical mechanics simulations for liquid hydrogen sulfide, methanethiol, ethanethiol, dimethyl sulfide, ethyl methyl sulfide, diethyl sulfide, and dimethyl disulfide. This also allowed detailed analyses of the structures for these liquids and study of condensed-phaseeffects for the conformational equilibria. The hydrogen bonding in the liquid thiols has been characterized, and the structures of liquid hydrogen sulfide and water have been contrasted.

Introduction There has been rapid progress in the application of computer simulation techniques to increasingly complicated and interesting chemical and biochemical systems.l A critical component of these molecular dynamics and Monte Carlo calculations is the potential functions that are needed to describe the interatomic interactions in the modeled systems. Great effort has gone into the development of intramolecular force fields for peptides, and several good alternatives are availablea2” The potential functions have largely been parameterized by using experimental vibrational data and structural results for crystals of amides and other organic molecules. The basis for the parametrization of peptide-water interactions has been limited, and it has been recognized repeatedly, in general, that “the most difficult set of parameters to derive a priori are the nonbonded To this end, we have been developing a complete set of potential functions for interpeptide and peptide-water interactions. The functions have been derived primarily by fitting directly to experimental thermodynamic and structural data on pure organic liquids and aqueous solution of small solutes representative of peptide backbones and side chains. This approach requires numerous fluid simulations; however, the resultant functions are clearly well-suited for modeling condensed-phase systems. The prior work on pure liquid water,’ hydrocarbons,s amides: and alcohols1° and on dilute aqueous solutions of amides,” hydrocarbons,12 carboxylate ions,13 and ammonium ions13 has yielded the parameters for 19 of the 24 most common peptide residues. The remaining ones are Cys, Met, His, Trp, and Arg. The work on sulfur compounds that covers the side chains for Cys, Met, and cystine bridges is presented here in the context of a general study of the properties and structure of liquid hydrogen sulfide, thiols, sulfides, and a disulfide. ~

(1) For reviews, see: (a) McCammon, J. A. Rep. Prog. Phys. 1984, 47, 1, (b) Computer Simulation of Chemical and Eiomolecular Systems; Beveridge, D. L.; Jorgensen, W. L., Eds. Ann. N.Y.Acad. Sci., in press. (2) For a review, see: Levitt, M. Annu. Rev. Eiophys. Eioeng. 1982, 11, 251. (3) Scheraga, H. A. Biopolymers 1981,20,1877. Nemenoff, R. A.; Snir, J.; Scheraga, H. A. J . Phys. Chem. 1978,82, 2504. (4) Gelin, B. R.; Karplus, M. Biochemistry 1979, 18, 1256. ( 5 ) Hagler, A. T.; Huler, E.; Lifson, S. J. Am. Chem. Soc. 1974, Y6, 5319. Hagler, A. T.; Lifson, S. Ibid. 1974, 96, 5327. (6) Weiner, S. J.; Kollman, P. A.; Case, D. A.; Singh, U. C.; Ghio, C.; Alagona, G.; Profeta, S.;Weiner, P. J. Am. Chem. SOC.1984, 106, 765. (7) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926. Jorgensen, W.L.; Madura, J. D.Mol. Phys. 1985,56, 1381. (8) Jorgensen, W. L.; Madura, J. D.; Swenson, C. J. J . Am. Chem. SOC. 1984, 106, 6638. (9) Jorgensen, W. L.; Swenson, C. J. J . Am. Chem. SOC.1985,107, 569. (10) Jorgensen, W. L. J. Phys. Chem. 1986, 90, 1276. (1 1) Jorgensen, W. L.; Swenson, C. J. J . Am. Chem. SOC.1985,107, 1489. (12) Jorgensen, W. L.; Gao, J.; Ravimohan, C. J . Phys. Chem. 1985,89, 3470. (13) Jorgensen, W. L.; Gao, J. J . Phys. Chem. 1986, 90, 2174.

0022-3654/86/2090-6379$01.50/0

TABLE I: ComputationalDetails of the Monte Carlo Simulations for Liquid Sulfur Compounds no. of configurations X liquid

H2S CH3SH CHSCHpSH (CH3)2S CHICHZSCH3 (CH3CH2)2S CH3SSCH 3

106 equil averaging 2.0 1.0 1.0 0.5 1.0 2.0 2.0

2.0 2.0 2.0 2.0 2.0 2.0 2.0

r,,

Ar,

9.0 10.0 11.0 11.0 12.0 14.0 12.0

0.22 0.22 0.20 0.25 0.15 0.15 0.14

A

A

A9,

A@, AV,

deg deg

A3

22 22 20 25 15 15 14

200 240 350 400 400 500 350

20

15 15 14

These first fluid simulations for liquid sulfur compounds also provide insights into many issues including hydrogen bonding in thiols, the influence of medium effects on molecular conformation, and the structure of liquid H2S vs. water. Computational Details Monte Carlo Simulations. Statistical mechanics simulations were carried out for seven liquids (hydrogen sulfide, methanethiol, ethanethiol, dimethyl sulfide, ethyl methyl sulfide, diethyl sulfide, and dimethyl disulfide) using potential functions described in the next two sections. Standard procedures were employed including Metropolis sampling, periodic boundary conditions, and the isothermal-isobaric (NPT) ensemble.s-10 Each system consisted of 125 or 128 molecules in a cubic cell. The external pressure was fixed at 1 atm, and the temperature was 25 O C except for the simulations of HIS and CH,SH which were run at the liquids’ The boiling points for boiling points, -60.34 and 5.96 0C.14,15 the other liquids range from 35.0 OC for C2H5SHto 109.6 OC for CH3SSCH,.I6.” Additional details are summarized in Table I. The intermolecular interactions were spherically truncated at cutoff distances, rcrbased on the S-S distances. A correction was made to the intermolecular energy for the Lennard-Jones interactions neglected beyond the cutoff. It amounts to ca. 2% of the total energy and was computed in a standard way.8 A correction was not made for the neglected electrostatic interactions; it is expected to be very small based on results for liquid water.’ New configurations of the systems were generated by random selection of a molecule followed by random translations, total rotations, and any internal rotations. The ranges of the three types of motion are given by (14) Giauque, W. F.; Blue, R. W. J. Am. Chem. SOC.1936, 58, 831. (15) Russell, H., Jr.; Osborne, D. W.; Yost, D. M. J. Am. Chem. SOC. 1942, 64, 165. (16) Haines, W. E.; Helm, R. V.;Bailey, C. W.; Ball, J. S.J. Phys. Chem. 1954, 58, 270. (17) Haines, W. E.; Helm, R. V.;Cook, G. L.; Ball, J. S. J . Phys. Chem. 1956, 60, 549.

0 1986 American Chemical Society

6380 The Journal of Physical Chemistry, Vol. 90, No. 23, 1986 TABLE 11: Geometrical and OPLS Parameters for Sulfur Compounds

S-H S-C S-S

C-C

Standard Geometrical Parametersb LHSH 92.0 LCSC LCSH 96.0 LCCS LSSH 98.0 LCSS

1.34 1.82 2.04 1.53

OPLS Parameters q, e u, A

atom or group S (HS) S (RSH) S (RSR)

s (R2S2) H (HzS) H (RSH) CH, (CH,SH) CH, (CH3SR) CHI (RCHZSR) CH3 (CHSSSR) CH2 (RCH2SSR)

-0.47 -0.45 -0.47 -0.30 0.235 0.27 0.18 0.235 0.235 0.30 0.30

c,

99.0 114.4 103.0

kcal/mol

3.70 3.55 3.55 3.55 0.0 0.0

0.25 0.25 0.25 0.25 0.0 0.0

a

a

3.80 3.80 3.80 3.80

0.17 0.118 0.17 0.118

“Standard alkyl group values from ref 8. bBond lengths in angstroms, bond angles in degrees.

f h r , *AB, and fAO in Table I. Attempts to change the volume of the systems were made every 700 configurations within ranges of f A V , and all intermolecular distances were scaled accordingly. The ranges were chosen to yield acceptance ratios of ca. 0.4 for new configurations. Each simulation had an equilibration phase prior to averaging over an additional 2 X lo6 configurations. The starting configurations were adopted from earlier simulations of liquid hydrocarbonss and alcohols.10 Convergence of the energy and density were readily achieved within the equilibration periods of 5 X lo5 to 2 x IO6 configurations. The statistical uncertainties reported here for computed properties are f 1u as obtained from separate averages over batches of 2 X lo5 configurations. The calculations were executed on a Gould 3218750 computer in our laboratory. Intermolecular Potential Functions. The optimized potentials for liquid simulations (OPLS) that we have been developing use the simple Coulomb plus Lennard-Jones form given in eq 1. The on a on b

AE,b = i

c (qiqf2/rV 4- Ajj/rijI2- Cij/riF)

(1)

J

interaction energy between molecules a and b is represented as the sum of the interactions between sites on the two molecules. The sites are located on the nuclei, except CH, groups are treated as single units centered on carbon. Standard combining rules are used for the short-range interactions such that Aij = (AijAjj)’/’ and Cij = (CiiC,,)1/2. The A and C parameters may also be expressed in terms of Lennard-Jones u’s and e’s as Aii = 4ciui12 and C, = 4eiu?. Standard bond lengths and bond angles are used for the monomers based on microwave data, as summarized in Table II.ls These are kept fixed during the simulations, though internal rotations are included (vide infra). The charges and Lennard-Jones parameters for the sulfur compounds were derived in our usual way for polar molecules by demanding reasonable results for bimolecular and ion-molecule complexes in the gas phase, consistency with the parametrization for other compounds, and successful description of the thermodynamics and structures of pure liquids.’-’’ ,The number of new parameters for the sulfur compounds is small since we insist on using the previously obtained Lennard-Jones parameters for CH, g r o ~ p s This . ~ ~ leaves ~ the Lennard-Jones parameters for sulfur and the charge distributions for the sulfur compounds to be determined. Though models were considered with interaction sites located in lone-pair positions on sulfur, no results were obtained that clearly demonstrated benefits of the added complexity. Thus, the simple OPLS model has an interaction site centered only on (18) Allinger, N. L.; Hickey, M. J. J . Am. Chem. SOC.1975, 97, 5167.

Jorgensen the sulfur nucleus. For reference, it may be noted that lone pairs are included on sulfur in the AMBER force field, while they are not in MM2.6J8J9 The final values for the parameters are summarized in Table I1 and are justified primarily in the Results section. However, a few comments on general patterns and comparisons with other potential functions are appropriate here. The Lennard-Jones u for sulfur is 3.55 A except in H2S for which a somewhat larger radius (3.70 A) was needed to obtain a reasonable density for the liquid. The former number is similar to the u of 3.56 8, used in AMBER, though the lone pairs on sulfur also have u’s of 1.07 A.6 The u for sulfur in the Gelin-Karplus potentials is smaller at 3.39 A and continues a pattern noted previously8 that would lead to overly dense fluids by use of these f ~ n c t i o n s .The ~ e for sulfur is difficult to optimize unambiguously and was fixed in the present work a t 0.25 kcal/mol. The value in AMBER is 0.20 kcal/mol, while the low estimate of 0.02 kcal/mol in the GelinKarplus potentials appears to come from underestimation of the polarizability of sulfur. Limited justification of the present value can be provided as follows. Optimization of Lennard-Jones parameters from several sources for noble gases indicates the e values increase by ca. 40% in going down a row in the periodic table.20 The e’s in the OPLS parameters for oxygens in water, alcohols, and amides are 0.15-0.21 kcal/mol, so an increment to 0.25 for sulfur seemed reasonable. The Lennard-Jones parameters for the CH, groups in the sulfur compounds were taken from the earlier work on hydrocarbonss and amide^.^ For the thiols as for alcohols,I0 the standard hydrocarbon parameters are used.8 However, for the sulfides and disulfides, somewhat smaller u’s are needed for the CH, groups adjacent to sulfur while more remote groups still use the hydrocarbon values.s Specifically, a u of 3.80 A is used for the adjacent groups, which is the same as was found desirable for N-alkyl groups in amidesg and recently for CH, groups cy to oxygen in ethers.21 The Lennard-Jones parameters for hydrogens on heteroatoms are taken to be zero in the OPLS format and in other potential functions for water and protein^.^,' Again, following the pattern for alcohols and amides?J0 charges have only been placed on sulfur and the two adjacent atoms. In the past, we had found that good initial estimates of charge distributions could be obtained from Mulliken populations for wave functions obtained by ab initio molecular orbital calculations with the 6-31G(d) basis set. However, this was not the case for sulfur; for example, the 6-31G(d) charge on sulfur in dimethyl sulfide is +O. 11. Thus, the charges were the principal parameters adjusted to obtain reasonable results for the bimolecular and ion-molecule complexes and pure liquids, as described below. The present model is so simple, however, that it should be realized that there is only one independent charge variable for H2S, sulfides, and disulfides and two variables for thiols. Dipole moments were also checked for thiols and sulfides; the experimental values are 1.5-1.6 D,22 while the OPLS charge distributions yield 2.2 and 2.7 D, respectively. An exact correspondence between experiment and the point charge models is not expected. For oxygen compounds, the OPLS charge distributions typically yield dipole moments 25% greater than experimental values. The larger differences for the sulfur compounds are undoubtedly somewhat associated with the greater length of CS than CO or CC bonds (Table XI). In view of the underestimate of heats of vaporization of sulfides (vide infra), even greater polarization seemed desirable; however, sodium ion affinities then become too large, so a compromise was necessary. For H2Sand dimethyl disulfide the point charge models yield dipole moments of 2.1 and 3.8 D, while the experimental values are 1.0 and 2.0 D.22*23 (19) Allinger, N. L.; Hickey, M. J.; Kao, J. J . Am. Chem. SOC.1976, 98, 2741. (20) Verlet, L.; Weis, J.-J. Mol. Phys. 1972, 24, 1013. (21) Jorgensen, W. L.; Briggs, J. M., to be submitted for publication. (22) CRC Handbook of Chemistry and Physics; 65th ed.; CRC: Boca Raton, FL, 1984; p E58. (23) Sutter, D.; Driezler, H.; Rudolph, H. D. Z . Nuturforsch., A: Astrophys., Phys. Chem. 1965, ZOa, 1676.

Potential Functions for Liquid Sulfur Compounds

The Journal of Physical Chemistry, Vol. 90, No. 23, 1986 6381 4 .o

TABLE III: Fourier Coefficients for Intramolecular Rotational Potential Functions"

ROTATIONAL ENERGY FUNCTION

~

molecule

bond

CH3CHZSH CH3CHZSCH3 (CHjCH2)ZS CHSSSCH,

C-S CH2-S C-S

S-S

Vo VI v2 V, 0.2100 -0.2067 -0.0733 1.2167 0.0000 1.4600 -1.0730 2.3100 0.0000 1.4600 -1.0730 2.3100 7.0100 1.5820 -8.6700 2.0080

" Vs in kcal/mol. 2 .oo

DIHEDRRL RNGLE DISTRIBUTIONS

FUNFTF]

ROTRTIONRL ENERGY

0

60

120

180

2r)O

300

IDERL GRS

1.2

360

DIHEDRAL RNGLE DISTRIBUTIONS .7En

-LIQUID - - _IDERL GRS

ETSH 25 C

0

60

120

180

240

300

360

PHI (DEG.1

Figure 2. Rotational energy function and dihedral angle distributions for ethyl methyl sulfide. Units as in Figure 1.

ROTATIONAL ENERGY FUNCTION 2

6.0

> 1

0

60

120

180

1

2L)O

I

300

360

PHI (DEG.1

Figure 1. Rotational energy function and dihedral angle distributions for ethanethiol. Units are kcal/mol for V(@)and mol %/deg for S(@).

.O 0

3 .oa

60

120

180

2110

300

360

DIHEDRRL RNGLE DISTRIBUTIONS --- LIQUID 2s c IDEAL GRS

ME-SS-ME

Internal Rotation. Torsional motion about the central C S bonds in CH3CHZSH,CH3CH2SCH3,and (CH3CH&S and about the SS bond in CH3SSCH3was included in the simulations. A review is available covering the implementation of the internal rotations and results for other The Fourier series in eq 2 is used V ( @ )= V, + 1/2V](1 + cos @)

2.25-

+ f/2V2(1 - cos 2@) + 1/2V3(1 + cos 3@) (2)

to represent the potential energy for rotation about a single bond. For diethyl sulfide, a separate Fourier series is used for each dihedral angle about the CS bonds and is augmented by a Lennard-Jones term between the methyl groups. As in the work on hydrocarbons and alcohols,sJO the Fourier coefficients were fit to results of M M 2 calculations for the rotational energy surf a c e ~ . ' ~This . ~ ~is well-justified by the extensive studies of sulfur compounds by Allinger and co-workers, who have demonstrated the accord with available experimental data.1s-'9 The Fourier coefficients are summarized in Table 111, and the rotational potentials for CH3CH2SH, CH3CH2SCH3,and CH3SSCH3 are illustrated in the top parts of Figures 1-3. The Lennard-Jones parameters for the 1-5 interaction in diethyl sulfide were taken from the work on hydrocarbons and are a = 4.0 A and = 0.0074 kcal/mol.* Ethanethiol has one trans (t) and two mirror-image gauche (g) conformers (Figure 1). From the potential function and the MM2 results,18trans is 0.21 kcal/mol higher in energy than gauche and the gauche-trans and gauchegauche barrier heights are 1.32 and 1.22 kcal/mol. IR and microwave studies concur that gauche is lower in energy by 0.2-0.3 kcal/mol, while the microwave study found g-t and g-g barriers of 1.41 and 1.76 k ~ a l / m o l . * ~ Ethyl methyl sulfide has a qualitatively similar profile (Figure 2), though now trans is lower than gauche. We fit to the MM2 relative energies for trans (O.O), the trans-gauche barrier (1.87), (24) Jorgensen, W. L. J . Phys. Chem. 1983, 87, 5304. (25) (a) Inagaki, F.; Harada, I.; Shimanouchi, T. J. Mol. Specrrosc. 1973, 46, 381. (b) Dung, J. R.;Bucy, W. E.; Wurrey, C. J.; Carreira, L. A. J. Phys. Chem. 1975, 79, 988.

0

60

120

180

290

300

360

PHI (OEG.1

Figure 3. Rotational energy function and dihedral angle distributions for dimethyl disulfide. Units as in Figure 1.

and the cis maximum (3.77). This results in the gauche form being 0.18 kcal/mol above trans, which is somewhat lower than the MM2 value of 0.29. As discussed elsewhere, the experimental preference has oscillated and the best available a b initio calculations indicate the two conformers are isoenergetic.26 For diethyl sulfide, the same Fourier coefficients were used as for ethyl methyl sulfide. The relative energies of the tt, tg, g+g+, and g+g- minima are 0.0, 0.18, 0.37, and 0.76 kcal/mol. The corresponding values from MM2 are 0.0, 0.35, 0.67, and 1.26 kcal/mol. The lower values with the present potential functions follow mostly from the gauche-trans discrepancy for ethyl methyl sulfide. The direction of the difference appears desirable based on the latest analyses for ethyl methyl sulfide.26 Spectroscopic evidence for at least three conformational forms of diethyl sulfide has been found in the gas and liquid states.27 Dimethyl disulfide only has two mirror-image gauche minima (Figure 3). The four Fourier coefficients were chosen to reproduce the MM2 energies of 0.0,7.01, and 10.60 kcal/mol for the gauche, trans, and cis rotamers and the location of the gauche minimum at a dihedral angle of 83.2O.I9 This is consistent with the available (26) Ohsaku, M.; Imamura, A. Mol. Phys. 1985, 55, 331. (27) Ohsaku, M.; Shiro, Y . ;Murata, H. B u f f .Chem. SOC.Jpn. 1973, 46, 1399.

6382 The Journal of Physical Chemistry, Vol. 90, No. 23, 1986 TABLE IV: Calculated Complexation Energies (kcal/mol) and Geometries complex linear (H2S), cyclic (H2S), bifurcated (H,S), antiparallel (H2S)2 linear (CH3SCH3)2 cyclic (CH,SCH3), bifurcated (CH,SCH,), antiparallel (CH,SCH3)2 linear (CH3SH), cyclic (CH3SH), linear (C2H5SH), HSH.*.OHZ HOH4H2 CH,SH*.*OH, HOH*.*SHCHI HOH**.S(CH,), Na+-.SH2 Na+-.SHCH, Na*.-S(CH3),

xy

SS SS

SS SS SS SS SS SS SS SS SS

SO OS SO OS OS NaS NaS NaS

Try, A ARCa 3.61 SSb 3.65 SSH 3.72 3.82 SSb 5.53 SSb 4.62 SSC 4.71 4.07 SSb 3.38 SSb 3.40 SSC 3.37 SSb 3.10 S o b 3.22 OSb 2.98 S o b 3.12 O S b 3.08 OSb 2.53 N a S H 2.47 N a S H 2.44 N a S C

eABC,

deg 150 40 63 126 55 71 147 40 146 139 150 142 140 121 134 142 130

-AEb 2.48 1.95 2.10 1.24 1.65 2.04 2.37 2.64 3.81 2.82 4.06 4.59 3.21 5.96 3.72 4.56 17.18 18.95 22.48

Jorgensen The linear form is also found to be favored for methanethiol dimer. The interaction energy is now -3.81 kcal/mol, and the SS separation is 3.38 8,. The hydrogen bond strengthens further for ethanethiol dimer (-4.06 kcal/mol), though the SS distance is essentially unaltered. These optimal distances do not appear so unreasonable when compared to some experimental data for condensed-phase systems. For example, X-ray structures for crystalline L-cysteine reveal SS contacts at distances as short as 3.68 Hydrogen bonding was also studied for linear complexes with water. HIS and CH3SH are both predicted to be better hydrogen-bond donors than acceptors. The interaction energies are -3.2 to -6.0 kcal/mol, and the SO distances are 3.0-3.2 A. These separations again appear desirable due to the observed SO disThe optimized value tances of ca. 3.4 A in crystalline ~-cysteine.~O of 3.08 h; for the SO separation in the water-dimethyl sulfide complex can also be compared with an SO distance of 3.0 A observed between the hydroxyl oxygen and sulfur for Tyr67 and Met80 in tuna ferrocytochrome c . ~ ' The greater hydrogen-bond donor than acceptor ability for the thiols vs. water is consistent with the formers' greater gas-phase acidities.32 We also performed ab initio calculations for 1 and 2 using the 6-31G(d) basis set with

a b is a point on the bisector of the RSR' or H O H angle of the hydrogen-bond acceptor. b-AE is the dissociation energy for the complex, M.-N M N.

-

+

experimental data including electron diffraction and microwave results of 83.9 f 0.9' and 84.7O for the dihedral angle.19-23,28 Finally, it should be noted that to enhance conformational transitions in the Monte Carlo simulations, we often use umbrella sampling methods." The barriers between conformers for ethanethiol and the sulfides are sufficiently low that this is not necessary at 25 OC. In view of the symmetry of the rotational potential for dimethyl disulfide, umbrella sampling was not used in this case as well.

Results and Discussion Bimolecular and Ion-Molecule Complexes. Geometry optimizations were carried out for various complexes of H2S, CH3SH, CH3CH2SH, and (CH3)2S,as summarized in Table IV. The linear, cyclic, bifurcated, and antiparallel forms for H2S dimers are the following: H Of\

H

/

/"\--s

H

H

linear

H

CYCllC

/

H

H

1

2

the monomers fixed in their experimental g e ~ m e t r i e s . ~The ~ computed interaction energy, r(SO), and 8 for 1 are -2.62, 3.63 A, and 153O, while they are -1.95, 3.81 A, and 108O for 2. Thus, the predicted order with the OPLS potentials appears correct, and the hydrogen bonds are stronger and shorter than gas-phase expectations. The antiparallel and bifurcated forms are predicted to be more favorable than the cyclic and linear forms for dimethyl sulfide dimer (Table IV). The lack of hydrogen bonding leads to greater interniolecular separations, and concomitantly, pure dipole-dipole interactions become more dominant. The optimal interaction energy is only -2.64 kcal/mol. A final check on the parameters was made by studying complexes with Na'. The Lennard-Jones parameters for Na+ were reported previously and used here with a charge of 1.0.34 Ab initio calculations for the three complexes in Table IV have been reported previously using the 6-31G(d) basis set with 3-21G optimized g e o m e t r i e ~ . ~We ~ also performed an MP2/6-3 1G(d) calculation for Na'SH, using 6-3 I G(d) optimized geometries; the complex has C,, symmetry with an N a S distance of 2.83 A and interaction energ of -14.5 kcal/mol. The corresponding OPLS results of 2.53 and -17.2 kcal/mol continue the pattern of somewhat accentuated interactions. For the CH3SH and (CH3)2Scomplexes, our best estimates of the interaction energies based on the ab initio results are ca. -17 and -20 kcal/mol with little shortening of the NaS separation. The OPLS predictions are ca. 2 kcal/mol too bound. Due to this observation, we were reluctant to increase the dipole moment of the sulfides further even though the heats of vaporization of the liquids are 7-15% too low as discussed in the next section.

+

x

In each case, the monomer geometries were kept fixed and typically one distance and the angle 6 were optimized by using the OPLS parameters. For H I S dimer, the linear form is found to be the most favorable with an interaction energy of -2.48 kcal/mol. The bifurcated dimer is less bound by 0.38 kcal/mol. The order agrees with the high level ab initio calculations including electron correlation Their predicted interaction (MP4/6-3 1+G(d,p)) of Frisch et energies are only -1.4 and -1.1 kcal/mol, and their SS separations of ca. 4.6 8, are significantly greater than our 3.6-3.7 8,. The need for shorter and stronger hydrogen bonds is typical of two-body potential functions for use in fluid simulations.' The lack of induced polarization in the potential functions may be particularly serious for sulfur. (28) Beagley, 8.; McAlloon, K. T. Tram. Faraday Soc. 1971,67, 3216. (29) Frisch, M. J.; Pople, J. A.; Del Bene, J. E. J . Phys. Chem. 1985,89, 3664.

(30) Harding, M. M.; Long, H. A. Acta Crystallogr., Sect B: Struct. Sci. 1968, 24, 1096. Kerr, K. A.; Ashmore, J. P.; Koetzle, T. F. Ibid. 1975, 31,

?n?? LULL.

(31) Takano, T.; Dickerson, R. E. J . Mol. Bioi. 1981, 153, 79. (32) Bartmess, J. E.; Scott, J. A.; McIver, R. T., Jr. J . Am. Chem. SOC. 1979, IOI, 6046. (33) (a) Binkley, J. S.; Whiteside, R. A,; Raghavachari, K.; Seeger, R.; DeFrees, D. J.; Schlegel, H.B.; Frisch, M. J.; Pople, J. A.; Kahn, L. R. "GAUSSIAN 82 Release H"; Carnegie-Mellon University: Pittsburgh, PA, 1982. (b) Francl, M. M.; Pietro, W. J.; Hehre, W. J.; Binkley, J. S.; Gordon, M. S.;DeFrees, D. J.; Pople, J. A. J . Chem. Phys. 1983, 77, 3054. (34) Chandrasekhar, J.; Spellmeyer, D. C.; Jorgensen, W. L. J. Am. Chem. Soc. 1984. 106. 903. (35) Smith,'S. F.; Chandrasekhar, J.; Jorgensen, W. L. J . Phys. Chem. 1982, 86, 3308.

The Journal of Physical Chemistry, Vol. 90, No. 23, 1986 6383

Potential Functions for Liquid Sulfur Compounds TABLE V Energetic Results (kcal/mol) for Liquid Sulfur Compounds

AH".. liquid H2S CHpSH CHpCHzSH (CH3)ZS CH3CH2SCHp (CHpCHM CH3SSCHp

-Ei

T, "C -60.34 5.96 25.0 25.0 25.0 25.0 25.0

Einm(B)

3.94 f 0.01 5.27 i 0.02 5.91 f 0.01 5.08 f 0.02 6.24 f 0.02 7.38 f 0.03 8.75 f 0.02

Eintm(1)

0.0 0.0

0.0 0.0

0.388

0.379 f 0.001

0.0

0.0

0.456 0.923 0.313

0.462 i 0.002 0.914 f 0.005 0.370 i 0.001

calcd

exptl

4.36 f 0.01 5.83 f 0.02 6.50 f 0.01 5.67 f 0.02 6.83 f 0.02 7.98 f 0.03 9.28 f 0.02

4.46' 5.87b 6.W 6.6Y 7.61d 8.55' 9.18c

'Reference 14. Reference 15. cReference 36. dReference 37. 'Reference 38. TABLE VI: Molecular Volumes and Densities for Liquid Sulfur Compounds

v,R3

T, OC -60.34 5.96 25.0 25.0 25.0 25.0 25.0

liouid --=---

HS CH3SH CH3CH2SH (CH3)2S CHpCH2SCH3 (CHpCHd2S CHpSSCHp

calcd 59.7 i 0.3 92.2 f 0.2 126.2 i 0.5 133.4 i 0.4 158.9 f 0.5 186.2 f 0.5 151.7 f 0.4

4 g/cm3 exDtl 58.51" 89.97b 123.8lC 122.47d 151.1Oe 180.14c 147.96c

calcd

exDtl

0.948 f 0.005 0.866 i 0.002 0.817 i 0.003 0.773 f 0.002 0.796 f 0.003 0.804 f 0.002 1.031 i 0.003

0.967" 0.888b 0.8332c 0.8423d 0.8368c 0.8312c 1 .0569c

'Reference 39. bReference40. CReference16. dReference 17. TABLE VII: Heat Capacities (cal/(mol deg)) and Isothermal Compressibilities (atm-I) for Liquid Sulfur Compounds"

T, liquid OC H2S -60.34 CHpSH 5.96 CH3CH2SH 25.0 (CH3)2S 25.0 CH3CH2SCH3 25.0 (CHpCH2)2S 25.0 CH3SSCHp 25.0

" Reference 4 1 .

Cpo(g) 8.0W 11.7Oe 17.37' 17.7lC 22.72f 27.95' 22.54'

CJ1) calcd exptl 18.0 i 25.1 f 28.9 f 28.6 f 34.9 f 39.1 f 31.4 f

1.1 1.5 1.5 1.3 1.6 1.4 1.3

1 0 6 ~calcd , 16.3b 227 i 30 21.28d 186 i 21 28.17' 176 f 24 28.23' 214 f 23 34.57f 192 f 26 40.97' 146 f 17 34.92e 81 f 9

Reference 14. Reference 42. Reference 36. /Reference 37. 'Reference 38.

Reference 15.

Thermodynamics. The thermodynamic results from the seven liquid simulations are compared with experimental data in Tables V-VII. The total energy of the liquid is given by the sum of the intermolecular, Ei(l), and intramolecular, Eintm(l),contributions. The relationship to the heat of vaporization is shown in eq 3 where M v a p = Eintra(g) - (Ei(1) + Eintn(1)) + R T (3) the intramolecular energy for the gas phase, Eintra(g),is obtained from a Boltzmann distribution over the internal rotational potential function, if present. The heat capacity, C,, for a liquid is estimated from the fluctuation in the intermolecular energy plus an intramolecular term taken as C, for the ideal gas less R. The isothermal compressibility is calculated from the volume fluctuations and has a substantial statistical uncertainty.' However, the total volume converges readily in Monte Carlo simulations in the NPT ensemble and is easily translated into a molecular volume and density. In Table V, the errors for the computed heats of vaporization are 1% or less except for the dialkyl sulfides. The error for dimethyl sulfide is 14.7%and lessens with increasing size to 10.2% for ethyl methyl sulfide and 6.7% for diethyl sulfide. I t may be (36) Wagman, D. D.; Evans, W. H.; Parker, V. B.; Schumm, R. H.; Halow, I.; Bailey, S. M.;Chumey, K. L.; Nuttall, R. L. J. Phys. Chem. R e j Data Suppl. 1982, lZ(2). (37) Scott, D. W.; Finke, H.L.;McCullough, J. P.; Gross, M.E.; Williamson, K. D.; Waddington, G.; Huffmann, H. M.J . Am. Chem. Soc. 1951, 73, 261.

(38) Scott, D. W.; Finke, H. L.; Hubbard, W. N.; McCullough, J. P.; Oliver, G. D.; Gross, M.E.;Katz, C.; Williamson, K. D.; Waddington, G.; Huffmann, H. M.J . Am. Chem. SOC.1952, 74, 4656. (39) Robinson, P. L.; Scott, W. E. J . Chem. Sac. 1936, 831. (40) Berthoud, A.; Brum, R.J. Chim. Phys. Phys.-Chim. Biol. 1924,22,

- ._. 1d l

(41) JANAF Thermochemical Tables, 2nd ed.;US.National Bureau of Standards: Washington, DC, 1971. (42) Barrow, G. M.; Pitzer, K. S.Ind. Eng. Chem. 1949,41, 2731.

TABLE VIII: Calculated Conformer Populations at compound conformer % gas' CHpCHzSH t 26.6 g 73.4 CH3CH2SCH3 t 43.1 g 56.9 (CH3CH2)2S t 47.5 52.5 B

25 OC % liquidb

27.4 f 0.6 72.6 f 0.6 42.5 f 1.0 57.5 f 1.0 48.5 f 0.7 51.5 f 0.7

" From Boltzmann distributions for the rotational potential functions. From the Monte Carlo simulations. noted that the largest errors in our study of liquid hydrocarbons were for propane.* As for dimethyl sulfide, an overly high volume is accompanied by a low heat of vaporization. Apparently, for the present model somewhat smaller Lennard-Jones u's are needed for the terminal methyl groups in these three-site systems than for larger molecules. Some improvement could also be obtained by using a smaller u and larger c for sulfur in sulfides than for thiols and disulfides or larger charge separation between sulfur and carbon. On balance, this introduction of additional parameters and potential worsening of the results for ion-dipole complexes, as mentioned above, seemed undesirable. The errors for the computed volumes and densities in Table VI average 2% except for the dialkyl sulfides for which the error is 6%. The computed volumes are all a little high with the largest error occurring for dimethyl sulfide, 8.9%. On this basis, a smaller u for sulfur would be beneficial; however, this would aggravate the description of the ion-molecule complexes. As shown in Table VII, the computed heat capacities are in good accord with experimental data. However, it should be realized that the ideal gas heat capacities less R make substantial contributions to C,, for the liquids, especially for the larger molecules. The computed isothermal compressibilities are also given in Table VII, though no experimental data could be found. The fluctuation properties converge relatively slowly in the simulations and have substantial statistical uncertainties. The magnitudes and trends for K seem reasonable in comparison with data for hydrocarbons; e.g., K declines from 175 X 10" to 124 X 10" atm-', in going from hexane to nonane at 25 0C.43 The particularly low K for dimethyl disulfide is consistent with its remarkably high heat of vaporization, 9.2 kcal/mol (Table V). From experimental data for alcohols and alkanes, we find a crude linear relationship between K and AHvap such that K = 310 19Mv,, where the units are lod atm-' and kcal/mol. The ex(43) Blinowski, A.; Brostow, W. J. Chem. Thermodyn. 1975, 7, 787. Eduljee, H. E.; Newitt, D. M.;Weale, K. E. J . Chem. SOC.1951, 3086.

Jorgensen

6384 The Journal of Physical Chemistry, Vol. 90, No. 23, 1986 5

DIHEDRAL ANGLE DISTRIBUTIONS IDEAL

1.2

0

60

120

180 2'+0 PHI IOEG.1

300

S-S RADIAL DISTRIBUTION FUNCTIONS

GAS

360

Figure 4. Dihedral angle distributions for rotation about the two C-S bonds in diethyl sulfide. Units as in Figure 1.

perimental K for liquid water and OUT predicted value for dimethyl disulfide deviate substantially from this equation toward the notably incompressible side. Overall, the simple OPLS model for sulfur compounds is found to provide a good description of the thermodynamic properties of liquid sulfur compounds. The average deviations from experiment for H2S, the thiols, and dimethyl disulfide are similar to the prior results for hydrocarbons, amides, and alcohols, while the errors for dialkyl sulfides are somewhat greater.*-1° Conformational Equilibria. The computed trans and gauche populations for liquid ethanethiol, ethyl methyl sulfide, and diethyl sulfide are compared with the gas-phase values in Table VIII. For diethyl sulfide, the trans and gauche populations are given averaged over both dihedral angles. The condensed-phase environment is found to have no statistically significant effect on the conformational equilibria. This is also the case for hydrocarbon liquids,8*44but small shifts in the conformer populations were found for alcohols.1° The intermolecular interactions are not strong enough for these three sulfur compounds to lead to conformational effects. Furthermore, in the OPLS model the dipole moments for these molecules are invariant to changes in the dihedral angles. The full dihedral angle distributions in Figures 1, 2, and 4 also show essentially no influence of the fluid environment. There is only a hint of some avoidance of dihedral angles nearing the cis rotamer for ethanethiol in Figure 1, an effect that is more pronounced for l-alkanols.'O In contrast, the dipole moment for dimethyl disulfide varies strongly with the dihedral angle, so a drift toward the cis rotamer is expected in the liquid phase since it is the most p0lar.2~This is clearly apparent in the dihedral angle distributions reported in Figure 3; the average value for the dihedral angle in the liquid is shifted 4 O toward cis from the gasphase reference point. On the technical side, the thorough sampling of the conformational space in the liquids is apparent in the symmetry of the dihedral angle distributions for ethanethiol, ethyl methyl sulfide, and diethyl sulfide (Figures 1, 2, and 4). However, the high barriers for dimethyl disulfide inhibit interconversion of the mirror-image conformers. In this case, the distribution of 6 0 g+ and 68 g- molecules at the end of the equilibration phase was maintained during the averaging and provided the asymmetry for s(@) in Figure 3. Umbrella sampling could possibly have reduced the asymmetry, though it also reduces the precision for all other computed quantities. Radial Distribution Functions-Thiols. Analysis of liquid structure often features radial distribution functions, gxy(r),that represent the probability of occurrence of atoms of type y a distance r from atoms of type x normalized for the bulk density of y atoms. The results from the Monte Carlo simulations are considered first for H2S and the thiols and then for the dialkyl sulfides and disulfide. The existence of any hydrogen bonding should be revealed in the S-S, S-H, and H-H radial distribution functions (rdfs) shown in Figures 5-7. The S-S rdfs for liquid CH3SH and CHzCH2SH (44) Kanesaka, I.; Snyder, R. G.; Straws, H.L. J . Chem. Phys. 1985,84, 395.

O I I

J 2

3

4

5 R(SS1.

6

7

8

9

h Figure 5. S-S radial distribution functions for liquid H2S and the thiols. Distances are in angstroms throughout. Successive curves offset 1.25 units along the ordinate.

i

DISTRIBUTION FUNCTIONS

I n 1

2

ETSH

3

4

5 R(SH1.

6

7

8

9

h Figure 6. S-H, radial distribution functions for liquid H2S and the thiols. Successive curves offset 1 unit along the ordinate. L)

0

H-H RRDIRL DISTRIBUTION FUNCTIONS

i

1

1

2

3

4

5 RIHHI.

6

6

7

8

9

Figure 7. H,-H, radial distribution functions for liquid H,S and the thiols. Successive curves offset 1 unit along the ordinate.

display distinct first peaks centered at 3.5 8, with heights of 1.7 and 1.9, respectively (Figure 5). Integration to 4 8, reveals 2 . 0 and 1.6 nearest neighbors for CH3SH and CH3CH2SH. (Note that the lower peak height, but larger integral for CH3SH, results from the higher number density of S atoms (Ns/ V) for CH3SH than CH3CHzSH.) Hydrogen bonding is suggested and is confirmed by the strong first peaks for the S-H rdfs in Figure 6. Their integrals to the minima near 3.05 8,are 0.9 for CH3SH and 0.8 for CH3CH2SH. Thus, it appears that the thiols have an average of a little less than one hydrogen bond per monomer and that there are contributions from non-hydrogen-bonded neighbors to the first peaks in the S-S rdfs. The same is true for the first peaks in the H-H rdfs (Figure 7) which integrate to 2.5 and 2.1 for CH3SH and CH3CH2SH out to 4.0 A. There is little structure beyond the first peaks in these rdfs. In comparison, liquid H2S has a strong first peak in the S-S rdf and relatively little structure in the S-H and H-H rdfs (Figures 5-7). This is consistent with substantial orientational disorder even at short distances and little firmly directed hydrogen bonding. The first peak in the S-S rdf integrates to 11.2 out to the minimum

The Journal of Physical Chemistry, Vol. 90, No. 23, 1986 6385

Potential Functions for Liquid Sulfur Compounds

r

ci-ci RRDIALDETRTGUTION FUNCTIONS 7

5

I ME-SS-HE

ET-S-ET

3 I

0

4

1

2

3

4

5 RlCCl.

6

7

8

9

A Figure 8. C , C , radial distribution functions for the liquid thiols. The curve for ethanethiol is offset 1 unit along the ordinate. 5

I

S-S RROIRL DISTRIBUTION FUNCTIONS

1

2

3

4

5 RISCI,

6

7

8

9

2

3

Y

5 R(CC1.

6

7

8

9

A Figure 10. S-methyl carbon radial distribution functions for the liquid sulfides and dimethyl disulfide. For ethyl methyl sulfide, the methyl carbon on sulfur is used. Successive curves offset 1 unit along the ordinate.

0 1

2

3

4

5

6

7

8

9

A Figure 9. S-S radial distribution functions for the liquid sulfides and dimethyl disulfide. Successive curves offset 1 unit along the ordinate. RISSI.

at 5.5 8,and is reminiscent of the rdfs for simple Lennard-Jones liquids. However, the separation of the first peak in the S-H rdf (Figure 6) is suggestive of some hydrogen bonding and integrates to 1.7 at 3.05 8,. Thus, the coordination number for sulfur is high as in structureless liquids, though there may be some hydrogen bonding embedded among many less orderly contacts with nearest neighbors. The remaining rdfs for the thiols provide few additional insights and are not shown here except for the Cl-Cl rdfs in Figure 8. In this case, the difference between CH3SH and CH3CH2SHis consistent with the occurrence of hydrogen-bonded chains or rings. The C1-0 rdfs for methanol and ethanol show the same pattem;1° the first band near 4 8, for CH3SH can be assigned to intermultimer contacts, while the bands near 5.3 8, for both thiols include the intramultimer contributions. As discussed below, the hydrogen-bonded multimers for the thiols are far smaller than for alcohols and, in fact, are mostly dimers and trimers. Radial Distribution Functions-Dialkyl Sulfides and Disulfides. Three rdfs for the dialkyl sulfides and dimethyl disulfide are shown in Figures 9-1 1. In general, there is little distinct structure in the rdfs for these liquids. The SS rdfs in Figure 9 have a broad first band that peaks near 5.2 8, in each case. A second band appears to be present centered near 6.5 8, for the disulfide, which likely contains contributions from the second, more distant sulfur of the neighboring monomers. Significantly shorter contacts are more common between the sulfurs and methyl groups (the methyl group on sulfur for CH3SCH2CH3),as shown in Figure 10. This results from normal packing effects and the favorable Coulombic interactions. The latter are actually less significant as witnessed by the similarity of the rdf for diethyl sulfide with the others in Figure 10 even though the methyl groups are uncharged in this case. The strongest first peak for the S-CH3rdfs is for dimethyl disulfide. The methyl groups of the disulfide have the largest positive charges and can be simultaneously proximate to a neighboring molecule's sulfur atoms due to the preferred gauche conformations. The rdfs between carbons a to the sulfur atoms are shown in Figure 1 1 and refer to the methyl carbons for CH3SCH2CH3. Packing effects are again apparent; the methylene groups for diethyl sulfide are clearly less accessible to each other than the terminal methyl groups for the other three systems.

0

I

4 1

d Figure 11. C , C , radial distribution functions for the liquid sulfides and dimethyl disulfide. Successive curves offset 1 unit along the ordinate. BClNDING ENERGY DISTRIBUTICINS

30

-25

-20

-15 -10 BONDING ENERGY

-5

0

Figure 12. Distributions of total intermolecular bonding energies for liquid H2S and the thiols. Bonding energies in kcal/mol. Units for the ordinate are mol %/(kcal/mol). In view of the lack of structure, integration limits for these rdfs are rather arbitrary. However, shallow minima are apparent near 4.5 8, in the S-CH3 rdfs in Figure 10. Integration to this point yields 3.4, 1.4, 2.1, and 3.2 neighboring methyl groups for the sulfurs in CH3SCH3,CH3SCH2CH3,(CH3CH2)2S,and CH3SSCH3. Since two methyl groups are identical by symmetry and can both contribute to the rdfs for these molecules except for CH3SCH2CH3,it appears that the interactions may be a little geometrically special with one or two neighbors for the dialkyl sulfides and with ca. three neighbors for the disulfide. Energy Distributions. The distributions in Figures 12-1 5 provide information on the energetic environment in the liquids. Figures 12 and 13 show the distributions of total intermolecular bonding energies for monomers in the liquids. The averages of these distributions are twice the Ei values in Table V. In each case, a unimodal profile is found covering a range of 10-15 kcal/mol for the bonding energies. The widths of the distributions are intermediate between those for alkanes* and alcohols10and

Jorgensen

6386 The Journal of Physical Chemistry, Vol. 90, No. 23, 1986

IBCINDING DISTRIBUTIONS 1

TABLE IX. Results of Hydrogen-Bond AIUI~YSCS for Liquid Thiols and H2S"

ENERGY

30

T, 'C no. of H bonds e(H bonds)

e(Cou1omb)

4-J) 4 deg

9, deg -25

-20

-15 -10 BONDING ENERGY

-5

5

CH$H

C2HSSH

5.96 1.01 -2.90 -2.76 -0.14 157 124

25.0 1.08 -3.07 -2.57 -0.49 155 122

7% of monomers in n H bonds

0

Figure 13. Distributions of total intermolecular bonding energies for liquid sulfides and dimethyl disulfide. Units as in Figure 12.

HS -60.34 0.08 -2.33 -2.48 0.15 118 130

0 1 2 3

92.0 7.8 0.2 0.0

27.0 47.0 24.0 2.0

24.4 45.6 27.4 2.6

"e's in kcal/mol. c(H bond) is the average hydrogen-bond energy which can be decomposed into the Coulomb, e(Coulomb), and Lennard-Jones, e(LJ), terms. A hydrogen bond is defined by an interaction energy of -2.25 kcal/mol or less. W v1

+

3

W V

0 ZI

d z

2 1

0

-2

-'I

0

2

DIMERIZRTION ENERGY

Figure 14. Distributions of individual interaction energies (kcal/mol) between molecules in liquid H2S and the thiols. Units for the ordinate are molecules/(kcal/mol). Successive curves offset 1 unit along the ordinate. 6

5

L) v)

W

-1

3 W V

d r

3

I

z

D

2

1

HE-S-ME 0

I

-6

-'I

!

-2

0

2

1

DIMERIZATION ENERGY

Figure 15. Distributions of individual interaction energies between molecules in the liquid sulfidts and dimethyl disulfide. Units as in Figure 14. Successive curves offset 1 unit along the ordinate.

reflect the strength and range of intermolecular interactions. When less variety is possible, the widths are narrower. Also, for the alcohols the distributions are bimodal due to distinction between molecules in only one hydrogen bond (chain ends) and the majority of molecules in two hydrogen bonds.I0 Such an effect is not apparent for H2S and the thiols in Figure 12. Energetic distinctions between molecules in different numbers of hydrogen bonds are not large in these cases, and the possible separate distributions merge together to yield the unimodal profiles. The energy pair distributions in Figures 14 and 15 show the distributions of individual molecule-molecule interactions in the liquids. The profiles for methanethiol and ethanethiol in Figure

14 have the classic bimodal form for hydrogen-bonded liquids. The bands at low energy correspond to the hydrogen-bonded neighbors, while the spikes between -1.5 and 1.O kcal/mol represent the weaker interactions with the many more distant molecules. Estimates of the number of hydrogen bonds per molecule can be obtained by integrating the low-energy bands to the inflection points near -2.25 kcal/mol, which is the same limit as used previously for analyses of the hydrogen bonding in liquid water.7 The results are 1.O and 1.1 hydrogen bonds for CH3SH and CH,CH,SH, which are a little higher than the estimates from the first peak in the S-H radial distribution functions (Figure 6). For comparison, the numbers of hydrogen bonds for water and methanol at 25 "C with the same cutoff are 3.6 and 2.0.7910 The results in Figure 15 for the dialkyl sulfides and dimethyl tiisulfide are typical for polar, aprotic liquids. There are some particularly favorable interactions with near neighbors extending to -4 kcal/mol for the disulfide. However, these molecules are not an energetically distinct group and merge into the continuum of weaker interactions. Clearly, these distributions spread with increasing size along the dialkyl sulfide series and with increasing polarity as in dimethyl disulfide. Both effects provide the opportunity for more favorable intermolecular interactions and lead to increasing heats of vaporization. Now, returning to Figure 14, the profile for H2S is an intermediate case. The clearly bimodal profile for hydrogen-bonded liquids is not present, though there is a shoulder at low energy. The optimal interactions of -2 to -2.5 kcal/mol are also sufficiently weak that designating them as hydrogen bonds is debatable. Integration to the -2.25 kcal/mol limit only yields 0.1 interaction. Overall, a possible description of liquid H2Sat its boiling point would be as a polar liquid with substantial orientational disorder and little hydrogen bonding. Hydrogen-Bonding Analyses. The hydrogen bonding in the liquid thiols and H2Swas further analyzed by using configurations saved every lo4steps during the Monte Carlo runs. A hydrogen bond was defined by the energetic cutoff of -2.25 kcal/mol described above. The results for the average numbers of hydrogen bonds, hydrogen-bond energies, and hydrogen-bond angles are recorded in Table IX. The hydrogen-bonding indices for the two thiols are similar. Both exhibit an average of one hydrogen bond per molecule and average hydrogen bond strengths of ca. 3 kcal/mol that are dominated by the Coulombic contributions. The average hydrogen bond angles O(S-H-S) and (P(H-43-H) are ca. 160° and 120°, which are close to the values for alcoholsI0 and are typical of a distribution of linear hydrogen-bonded geometries. The percentages of molcules in 0, 1, or 2 hydrogen bonds are roughly 25, 50, and 25. This is consistent with a large fraction of linear trimers since they require a 2: 1 ratio of molecules in one and two hydrogen bonds. Stereoplots of configurations from the simulation do show trimers as well as dimers and monomers. The substantial per-

Potential Functions for Liquid Sulfur Compounds

The Journal of Physical Chemistry, Vol. 90, No. 23, 1986 6387

Figure 16. Stereoplots of configurations from Monte Carlo simulations of liquid water at 25 “C (top) and H2S at -60 molecules in both plots.

O C

(bottom). There are 125

centage of monomer found here is notable in comparison to the violate the spirit of hydrogen bonding and are perhaps better 0-2% predicted for liquid alcohols.I0 Of course, the “monomers” described simply as dipolar interactions. For comparison, stermay have several interactions just above the hydrogen-bonding eoplots of configurations from the Monte Carlo simulations of cutoff since it is clear from Figure 12 that the “monomers” must liquid H2S at -60 OC and liquid TIP4P water’ at 25 OC are still be bound by 5-10 kcal/mol in the liquids. presented in Figure 16. In viewing the plots, the periodicity should The major conclusion here is that there is substantial hydrogen be recalled so that molecules near one face of the cube are bonding to the liquid thiols dominated by the Occurrence of linear proximate to molecules near the opposite face and the edges of dimers and trimers. This may be surprising in view of some the illustrated cube are somewhat outside the bounds of the pequalitative notions in the literature. For example, in the book by riodic cube. Another point for reference is that the actual volume Vinogradov and Linnell, Met is listed as a hydrophobic group that per molecule for H2S is twice that for H 2 0 though the stereoplots does not participate in hydrogen bonds, while Cys is possibly a have been drawn the same size for ease of viewing; correcting for hydrogen-bond donor, but probably not a hydrogen-bond accepthis would compress the plot for water and make it appear much more dense. For water, the hydrogen bonding is clear with at least t0r.4~ Also, Pimentel and McClellan note that the S-H group three readily identifiable hydrogen bonds for most molecules. shows hydrogen bonding with strong bases, e.g., amines; however, they state that “the relative weakness of S-H as a proton donor There are also no obviously repulsive neighboring pairs; the accounts for the absence of H bonds in some systems.”46 From majority of oritentations for neighbors are for bent hydrogen bonds the present results, it is predicted that the sulfur functionalities of the linear form. On the other hand, for H2S there is much orientational disorder. There are some neighboring pairs in hyin Cys and Met are perfectly capable of participating in hydrogen drogen-bonding orientations, but their Occurrence is much more bonds. This position is also supported by the crystal structures ~ ~ ~ than for water. The number of neighbors in the first for L-cysteine and tuna ferrocytochrome c noted a b o ~ e . ~ ~random Futhermore, ab initio calculations and our potential functions both solvent shell is significantly larger for HIS and not as easily identifiable as the 3-4 hydrogen-bonded neighbors for water. find the S-H group to be a stronger hydrogen-bond donor than acceptor (Table IV). Finally, I R studies of liquid thiols have been There is also a significant number of bifurcated interactions for H2S as illustrated by the five molecules at the bottom of the box. undertaken and provide clear evidence for extensive association In inert gas matrices, methanethiol due to hydrogen bonding!’ The bifurcated form is favored on a dipole-dipole basis and is only 0.38 kcal/mol less favorable than the linear form for H2S (Table and ethanethiol also show multimer formation even at very low IV), while the difference is ca. 2.5 kcal/mol for several potential concentrations; extensive hydrogen bonding is found as well in functions for the water dimer including TIP4P.’ On the experthe pure solids with band splittings due to the presence of both imental side, vibrational spectroscopic studies for liquid H2S do gauche and trans conformers for e t h a n e t h i ~ l . ~ ~ For liquid H2S at -60 “ C the situation is less clear. If one uses not seem to have been performed; however, there is IR evidence for H2S multimer formation in inert gas matrices at 20 K with the same energetic criterion as for the thiols, there is only 0.1 hydrogen bond per molecule and 92%monomer as listed in Table a preference for the linear dimer.49*50 IX. If the hydrogen-bond limit were relaxed to 2 or 1 kcal/mol, Conclusion then there would be 0.2 or 1.0 hydrogen bond per molecule. Our subjective opinion is that interactions weaker than 2 kcal/mol Intermolecular potential functions have been developed for use in computer simulations of thiols, sulfides, and disulfides including the sulfur-containing side chains of proteins. The good accord (45) Vinogradov,S. N.; Linnell, R. H. Hydrogen Bonding, Van Nwtrand with experimental data for the thermodynamic properties of the Reinhold: New York, 1911; p 234. (46) Pimentel, G. C.; McClellan, A. L. The Hydrogen Bond; W. H. Freeman: San Francisco, 1960; p 201. (47) Bicca de Alencastro, R.; Sandorfy, C. Can. J. Chem. 1972,50,3594. (48) Barnes, A. J.; Hallam, H. E.; Howells, J. D. R. J. Chem. SOC.,Far-

aday Tram. 2 1972,68, 131.

(49) Tursi, A. J.; Nixon, E. R. J . Chem. Phys. 1970, 53, 518. (50) Barnes, A. J.; Howells, J. D. R. J. Chem. Soc., Faraday Trans. 2 1912, 68, 129.

6388

J. Phys. Chem. 1986, 90, 6388-6392

liquid sulfur compounds and the reasonable description of bimolecular and ion-molecule complexes obtained with the OPLS potentials are impressive in view of the simple form and parametrization of the functions. Insights into the structures of the liquids were also obtained. Hydrogen bonding is evident for liquid methanethiol and ethanethiol; there is a substantial occurrence of linear dimers and trimers with an average of one hydrogen bond per molecule. The existence of hydrogen bonding in liquid H2S at its boiling point is questionable. However, it is reasonable to describe this liquid and the liquid sulfides and dimethyl disulfide as dipolar fluids with substantial orientational disorder. The effect of the condensed phase on the conformational equilibria for

ethanethiol, ethyl methyl sulfide, and diethyl sulfide was also studied and found to be negligible. For dimethyl disulfide, a small shift is predicted in the average dihedral angle favoring a more polar form in the liquid than in the gas phase.

Acknowledgment. Gratitude is expressed to the National Institutes of Health and the National Science Foundation for support. Assistance with plotting programs was provided by Jiali Gao. Registry No. H2S,7783-06-4;CH,SH, 74-93-1; CH3CH2SH,7508-1;(CH3)2S,75-18-3;CH,CH2SCH, 624-89-5;(CH,CH2)2S,35293-2;CHSSSCHS,624-92-0.

Surface-Enhanced Resonance Raman Spectra of Adriamycin, 11-Deoxycarminomycin, Their Model Chromophores, and Their Complexes with DNA Giulietta Smulevich* and Alessandro Feis Dipartimento di Chimica dell’ Universitri, 501 21 Firenze, Italy (Received: June 2, 1986)

The surface-enhand resonance Raman spectra on Ag sols of adriamycin, 1 1-deoxycarminomycin,their model chromophores 1,4-and 1,8-dihydroxyanthraquinones,and their complexes with DNA were measured. The well-detailed spectra allowed us to obtain, by the combined analysis in terms of symmetry and pseudosymmetry, a nearly complete vibrational assignment of the resonance Raman active modes. The spectral perturbations induced by the adsorption of the compounds onto the Ag particles, by comparison with their resonance Raman spectra in solution, were explained in terms of interaction between one C=O.-H-O group of the chromophore and the silver surface. The intensity reduction of some bands associated with the HOCCC==O groups observed in the drug/DNA complexes was interpreted in terms of changes between the ground and the excited states of the normal coordinates and/or their equilibrium positions. The inferred structures of the complexes were found to be consistent with intercalation between daunorubicin and the DNA fragment d(CpGpTpApCpG) already obtained from X-ray measurements.

Introduction Glycoside anthracycline antibiotics are widely used in the treatment of cancer.’ Following adriamycin (hereafter ADM) and daunorubicin (DNR), new series of glycoside anthracycline antibiotics have been developed in order to improve their efficacy and minimize the toxic effects.2 The biological activity of all these drugs has been attributed to the formation of the intercalation complexes between the chromophore framework and the base pairs of DNA.) Much research has been carried out in order to elucidate the interaction mechanism, in particular resonance Raman spectroscopy has been widely used because of its selectivity which permits the observation of only bands due to the vibrations of the chromophoric framework. Such studies have been reported for adriamycin?” the adriamycin-DNA complex,6.’ and copper- and iron-adriamycin complexes.8-10 A handicap in obtaining highly detailed resonance (1) Remers, W. A. The Chemistry of Antitumor Antibiotics; Wiley: New York, 1979; Vol. 1. (2) Arcamone, F. Med. Res. Reu. 1984,4, 153-188 and references therein. (3)Aubel-Sadron, G.; Londos-Gagliardi, D. Biochimie 1984,66, 333-352 and references therein. (4) Hillig, K. W.; Morris, M. D. Biochem. Biophys. Res. Commun. 1976, 71, 1228-1233. (5) Manfait, M.; Bcmard, L.; Theophanides, T. J . Raman Spectrosc. 1981, 1 1 , 68-74. (6) Angeloni, L.; Smulevich, G.; Marzocchi, M. P . Spectrochim. Acta, Part A 1982, 38A, 213-217. (7) Manfait, M.; Alix, A. J. P.; Jeannesson, P.; Jardillier, J.-C.; Theophanides, T. Nucleic Acids Res. 1982, 10, 3803-3816. (8) Beraldo, H.; Gamier-Suillerot, A,; Tosi, L.; Znorg. Chem. 1983, 22, 41 17-41 24. (9) Beraldo, H.; Gamier-Suillerot, A.; Tosi, L.; Lavelle, F. Biochemistry 1985, 24, 284-289.

0022-3654/86/2090-6388$01.50/0

Raman spectra is represented by the strong fluorescence of these compounds. Surface-enhanced Raman scattering (SERS) and surface-enhanced resonance Raman scattering (SERRS) have been used to obtain information on biomolecules.1’J2 In particular it is a powerful method for fluorescent chromophores. In fact it allows total quenching of the chromophore fluorescence and at the same time magnifies the normal Raman scattering by several orders of magnitude, giving rise to structural information not always obtained with other methods. In order to take advantage of the SERS effect it is necessary to understand the mechanism of binding of adsorbed molecules to metal surfaces. In particular the possibility must be considered of structural alteration induced by surface forces. It was demonstrated by recent studies on different biomolecules that while for flavin chromophores the enzyme active site remains intact,13 for heme proteins the SERS spectra, produced by heme protein adsorbed on silver colloids, are due to surface-bonded p-oxobridges iron porphyrin dimers, formed by disruption of the heme-binding pockets of the p r ~ t e i n . ’ ~ In the present study we examine the SERRS spectra on Ag colloids obtained for highly fluorescent antitumor drugs, their complexes with DNA, and their model chromophores. (IO) Dutta, P. K.; Hutt, J. A. Biochemistry 1986, 25, 691-69s. (1 1) Koglin, E.; Sequaris, J. M. in Spectroscopy of Biological Molecules; Alix, A. J., Bernard, L., Manfait, M., Eds.; Wiley: New York, 1985; pp 221-225 and references therein. (12) Cotton, T. M. in Surface and Interfacial Aspects of Biomedical Polymers; Andrade, J. D., Ed.; Plenum: New York, 1985; Vol. 11. (13) Copeland, R. A.; Fodor, S.P. A.; Spiro, T. G. J. Am. Chem. SOC. 1984, 106, 3872-3874. (14) Smulevich, G.; Spiro, T. G. J . Phys. Chem. 1985, 89, 5168-5173.

0 1986 American Chemical Society