Chlorinated ethanes in aqueous solution: parameterization based on thermodynamics of hydration

CHEMICAL

30 August 1996

PHYSICS LETTERS ELSEVIER

Chemical Physics Letters 259 (1996) 142-145

Chlorinated ethanes in aqueous solution: parameterization based on thermodynamics of hydration M.D. Paulsen, T.P. Straatsma Environmental Molecular Sciences Laboratory, Pacific Northwest National Laboratory, Richland, WA 99352, USA

Received 10 May 1996

Abstract

Lennard-Jones parameters for short chain alkanes and chloroethanes were optimized in order to reproduce experimental absolute hydration free energy data. Absolute solvation free energies were calculated using the recently developed separation-shifted scaling method in conjunction with thermodynamic integration calculations. Coefficients for the repulsive interaction of carbon, hydrogen, and chlorine with SPC/E water were optimized in order to simultaneously improve the agreement between calculated and experimental absolute hydration free energies for both unsubstituted short chain alkanes and chloroethanes. However, the optimized parameters resulted in substantially poorer agreement for pure liquid densities of chloroethanes.

1. Introduction One of the important methodological advancements in the field of molecular simulations has been the development of techniques for the computation of free energies. These techniques permit the calculation and comparison with experimental data of a variety of thermodynamic properties, such as the thermodynamics of solvation. To date, the vast majority of such studies have employed molecular mechanics force fields rather than quantum mechanical descriptions for both solutes and solvents. Several different force fields appropriate for organic or biochemical systems have been developed [1-5]. These force fields have been optimized to reproduce a variety of data including crystallographic and spectroscopic data, ab initio calculations, and properties of pure liquids such as densities and heats of vaporization. However, they have not generally been pa-

rameterized against free energies such as the free energy of hydration. In this work, free energy calculations have been used to evaluate the suitability and perform optimization of force field parameters to reproduce absolute free energies of hydration of small alkanes and chlorinated ethanes.

2. Methods Free energy differences can be evaluated from molecular dynamics simulations using the technique of thermodynamic integration [6-8]. The two states, A and B, for which the free energy difference is to be computed are coupled by a control variable A which is slowly varied during the course of the molecular dynamics trajectory. For hydration free energies, state A is a gas phase solute and pure solvent which do not interact with each other while

0009-2614/96/$12.00 Copyright © 1996 Elsevier Science B.V. All rights reserved PII S0009-261 4(96)00729-4

M.D. Paulsen, T.P. Straatsma / Chemical Physics Letters 259 (1996) 142-145

Table 1 Alkane free energyof solvation OPLS

Revised

ethane

(kJ/mol) 14.24±1.43

(kJ/mol) 7.38+1.26

Experimentala (IO/mol) 7.66

propane butane

15.11 + 1.66 17.44 4- 1.66

9.37 + 1.29 6.43 4- 1.64

8.18 8.70

a Ref. [15].

constant pressure of 1 × 10 s Pa. The pressure and temperature of the system were maintained by coupling to an external bath [13]. All bonds involving hydrogen were constrained using SHAKE [14], and a 2.0 fs time step was used in all simulations.

4. R e s u l t s a n d d i s c u s s i o n

state B is a solution. The free energy difference between states A and B where state A is given by A = 0 and state B by A = 1 is then given by AG = f

143

(OH(A)/0A) dA = ~ (0H(A)/0A)AA. i This expression requires that the derivative of the Hamiltonian with respect to the control variable A be evaluated at each time step and then an ensemble average is constructed for each discrete value of A used in the trajectory, "o

3. Computational details For the determination of hydration free energies, multiconfiguration thermodynamic integration in conjunction with the separation-shifted scaling method as implemented in the ARGOS simulation package was used [9-11]. For all solute simulations, 21 equally spaced values of A were used. At each value of A, 750 steps of equilibration and 600 steps of data gathering were performed. The model system was a periodic box, 2,4 nm on a side containing the solute and 410 solvent molecules represented using the S P C / E water model [12]. All simulations were performed at a constant temperature of 298 K and a

Initial calculations of the hydration free energy for hexachloroethane (8 _+ 2 kJ/mol) and 1,1,1-trichloroethane (5 _ 2 kJ/mol) using OPLS nonbonded parameters gave very poor results for both absolute and relative free energies. To further investigated the source of this discrepancy, multiconfiguration thermodynamic integration calculations were performed to evaluate absolute hydration free energies for unsubstituted alkanes using OPLS nonbonded parameters [16]. The results are summarized in Table 1. For the solutes examined, ethane, propane, and butane, the calculated hydration free energies are substantially too positive. However, since the error is systematic, calculated relative hydration free energies for the three solutes are in better agreement with experiment than the absolute hydration free energies. In order to optimize, the agreement between calculated and experimental free energies of hydration, the derivative of the energy with respect to each of the nonbonded force field parameters was determined. This analysis indicated, as expected for nonpolar solutes, that the free energies of hydration for these solutes are relatively insensitive to small changes in the partial charge assignment. However, as the derivatives summarized in Table 2 indicate, the calculated hydration free energies will be sensitive to small variations in the repulsive portion of the Lennard-Jones term of the force field. The deriva-

Table 2 Energy derivatives with respect to C~2 nonbonded parameters Ethane

Propane

Butane

(aU/~CI2)C_ o (nm- 12) (~U/aCi2)n_ o (nm- 12)

3.0240 x 106 5.5221 X 107

3.6069 × 106 6.6847 x 107

4.0548 X 106 7.4369 )< 107

(8U/~Ci2)c1_ o ( n m - t z )

HCA 9.3377 × 106

TetCA 6.3339 )< 106

144

M.D. Paulsen, T.P. Straatsma / Chemical Physics Letters 259 (1996) 142-145

Table 3 Chloroalkane free energyof solvation Initial Revised

Experimental a

(kJ/mol)

(kJ/mol)

(kJ/mol)

HCA

9.59 + 2.37

- 4.89_+ 1.97

TetCA

5.295:1.86

-8.725:1.81

- 5.88 -9.87

a Ref. [15].

tives in Table 2 were used in combination with the discrepancies between calculated and experimental free energies summarized in Table 1 to obtain revised r 12 coefficients for carbon and hydrogen. Since the derivatives summarized in Table 2, are strictly valid only at the values of the nonbonded parameters at which they were calculated, it was necessary to repeat the optimization procedure three times before obtaining satisfactory agreement between calculated and experimental hydration free energies. The calculated hydration free energies using the revised nonbonded parameters for carbon and hydrogen are shown in Table 1. The difference between calculated and experimental hydration free energies for ethane and propane is now less than the estimated uncertainty in the calculated values. The discrepancy between calculated and experimental values for butane is slightly larger but still is significantly smaller than with the original OPLS parameters. The revised hydrogen and carbon parameters were then used to evaluate hydration free energies for two chloroalkanes, hexachloroethane and 1,1,2,2-tetrachloroethane. As shown in Table 3, even with revised carbon and hydrogen parameters, there are still

large discrepancies between calculated and experimental solvation free energies for these two solutes. As in the case of the unsubstituted alkanes, the calculated hydration free energy for hexachloroethane and 1,1,2,2-tetrachloroethane will not be particularly sensitive to small variations in the partial charge assignments. Thus improved agreement with experiment was achieved by optimization of the r j2 coefficient for chlorine. After two rounds of optimization, the agreement between calculated and experimental hydration free energies was significantly improved as shown in Table 3. The difference between calculated and experimental values for both solutes is now of the same order as the estimated uncertainty in the calculated values. To further test the applicability of the revised nonbonded carbon, hydrogen, and chlorine parameters, hydration free energies for all possible chloroethanes were calculated. The results of these chloroethane calculations are presented in Table 4 and compared with the available experimental data. For chloroethane, 1,1-dichloroethane, 1,1,2-trichloroethane, and pentachloroethane, the use of the revised carbon, hydrogen, and chlorine resulted in excellent agreement between experimental and calculated hydration free energies. For the other two solutes, 1,2-dichloroethane and 1,1,1-trichloroethane, use of the revised Lennard-Jones parameters resulted in calculated hydration free energies which were too low; - 9 . 2 6 kJ/mol for 1,2-dichloroethane and - 5 . 2 2 kJ/mol for l,l,l-trichloroethane. Since both of these compounds possess significant dipole moments in

Table 4 Chloroethane free energies of hydration and pure liquid densities Solute

AG ealc ( k J / m o l )

AG exp ( k J / m o l ) a

MCA 1,2-DCA 1,1-DCA 1,1,2-TCA 1,1,I-TCA I,I,I,2-TetCA 1,1,2,2-TetCA PCA HCA

-1 5:2 -6 _ 2 -5 5:2 - 10 5 : 2 - 3 5:1 - 10 4- 2 - 9 5:2 - 7 5:2 -5 5:2

-2.5 -6.1 -3.5 - 8.4 - 1.1 - 9.9 - 5.7 -5.9

Pola ( k g / m - 3) 823 b 1110 1220 1380 1520 1550

a Data taken from Ref. [15]. b Simulations run at 273.15 K. All other pure liquid simulations were run at 298.15 K c Experimental data taken at 273.15 K. All other densities were determined at 298.15 K

Prey ( k g / m - 3)

Pexp ( k g / m - 3)

1100 b 1429 1430 1680 1666 1897 1840 2033 2221

924 c 1245 1176 1439 1339 1541 1600 1679 2091

M.D. Paulsen, T.P. Straatsma / Chemical Physics Letters 259 (1996) 142-145

145

dilute aqueous solution, improved agreement with experiment can be obtained by modest scaling of the partial charges. The hydration free energies for these two solutes, calculated using the scaled partial charges are included in Table 4.

results in much better agreement with experimentally determined thermodynamic data.

The initial force field parameters used for the optimization with absolute free energies of hydration as criterion were taken from a set optimized to reproduce pure liquid properties. After the optimization with emphasis on the thermodynamics of hydration, it can be expected that the optimized parameters will do less well for the description of pure liquid properties. To determine to what extent the optimized parameters are still able to reproduce pure liquid interactions, molecular dynamics simulations of 100 ps in length were carried out for each of the chlorinated ethanes. The parameters for the pure liquid simulation were obtained by using the

Pacific Northwest National Laboratory is a multiprogram national laboratory, operated for the U.S. Department of Energy by Battelle Memorial Institute under Contract DE-AC06-76RLO 1830. This work was supported in part by a grant (KP0402) from the Health Effects and Life Sciences Research Division, and by the Environmental Molecular Sciences Laboratory construction project, Office of Health and Environmental Research of the Office of Energy Research of the U.S. Department of Energy.

Lorentz-Berthelot combination rules. The density of the liquid, which had been used as one of the optimization criteria for the original parameters, are listed in Table 4 together with the experimentally determined values. Using the optimized parameters, the calculated densities are systematically too high.

[1] W.F. van Gunsteren and H.J.C. Berendsen, Groningen Molecular Simulations (GROMOS)Library Manual (Binmos, Groningen, 1987). [2] B.R. Brooks, R.E. Bmcoleri, B.D. Olafson, D.J. Slater, S.

5. Conclusion In this article, a force field reparameterization is presented for chlorinated ethanes in which the thermodynamics of hydration has been used as the criterion for optimization. The incentive for this study has been to develop parameters for the chlorinated ethanes for use in molecular simulation studies of these compounds in aqueous environment, in which the calculation of thermodynamic properties, in particular free energy differences, are of key importance. The analysis of parameterizations in this article illustrates that, at least, for effective pair models of molecular interactions, the simple LorentzBerthelot combination rules may be inadequate. The model parameters for small alkanes and chlorinated ethanes give satisfactory results for the pure liquid simulations. The interaction with water as a solvent is however, not properly described using the simple combination rules. The reparameterization performed

Acknowledgement

References

Swaminathan and M. Karplus, J. Comp. Chem. 4 (1983) 187.

[3] S.J. Weiner, P.A. Kollman, D.T. Nguyen and D.A. Case, J.

Comp. Chem. 7 (1986) 230. [4] W.L. Jorgensen and J. Tirado-Rives, J. Am. Chem. Soc. 110 (1988) 1657. [5] W.D. Comell, P. Cieplak, C.I. Bayly, I.R. Gould, K.M. Merz Jr,, D.M. Ferguson, D.C. Spellmeyer, T. Fox, J.W. Caldwell and P.A. Kollman, J. Am. Chem. Soc. 117 (1995) 5179. [6] W.L. Jorgensen, AccountsChem. Res. 22 (1989) 184. [7] T.P. Straatsma and J.A. McCammon, in: Methods in enzymology, ed. J.J. Langone (Academic Press, San Diego, 1991)

pp. 497-511. [8] P. Kollman,Chem.Rev. 93 (1993)2395. [9] T.P. Straatsma and J.A. McCammon, J. Comp. Chem. 11 (1990) 943.

[10] M. Zacharias,T.P. Straatsmaand J.A. McCammon, J. Chem. Phys. 100 (1994) 9025. [11] T.P. Straatsma and J.A. McCammon, J. Chem. Phys. 95

(1991) 1175. [12] H.J.C. Berendsen, J.R, Grigera and T.P. Straatsma, J. Phys. Chem. 91 (1987) 6269. [13] H.J.C. Berendsen, J.P.M. Postma, W.F. van Gunsteren, A. DiNola and J.R. Haak, J. Chem. Phys. 81 (1984) 3684. [14] J. Ryckaert, G. Ciccotti and H.J.C. Berendsen, J. Comput.

Phys. 23 (1977)327. [15] S. Cabani, P. Gianni, V. Mollica and L. Lepori, J. Solution

Chem. 10 (1981) 563. [16] W.L. Jorgensen, J.M. Briggs and M,L. Contreras, J. Phys.

Chem. 94 (1990) 1683.