12162
J. Phys. Chem. 1996, 100, 12162-12171
Molecular Dynamics Simulation of a Micellar System: 2,3,6,7,10,11-Hexakis(1,4,7-trioxaoctyl)triphenylene in Water Tim Bast* and Reinhard Hentschke Max-Planck-Institut fu¨ r Polymerforschung, Postfach 3148, 55021 Mainz, Germany ReceiVed: December 20, 1995; In Final Form: March 6, 1996X
We present the results of an atomistic molecular dynamics simulation based on the AMBER/OPLS force field applied to segments of isolated one-dimensional micelles, 2,3,6,7,10,11-hexakis(1,4,7-trioxaoctyl)triphenylene, in aqueous solution using the SPC/E water model. The quantities which we study include the intramicellar monomer structure and dynamics, e.g., the equilibrium monomer-monomer separation along the micelle as well as the rotational dynamics of the monomers. We also study the micelle-water interface, which yields the effective micellar diameter, and the flexibility of the micelle in terms of its persistence length as a function of temperature. In addition, we determine the micelle size distribution at low concentrations via the free enthalpy gain per monomer-monomer contact using a hydration shell model in combination with thermodynamic integration.
I. Introduction Many amphiphilic molecules in aqueous solution spontaneously self-assemble into labile anisotropic aggregates which exhibit complex liquid-crystalline phase behavior (e.g., ref 1). Over the past decade, theoretical modeling of the interplay between self-assembly and lyotropic phase behavior in these systems has been addressed by various authors (e.g., refs 2 (review article) and 3-6). What is common to these models is that they focus on rod-shaped micelles (some exceptions are discussed in ref 2) and that they are extensions of Onsager’s excluded volume theory for the isotropic-to-nematic transition in a dilute systems of slender rods. The theoretical models usually depend on a number of a priori unknown parameters such as the monomer shape and size, the flexibility of the micelle expressed in terms of its persistence length (for a rodlike micelle), the monomer-monomer contact enthalpy within a micelle, etc. These parameters are either adjusted by fitting the model results to the experiment or deduced from independent experimental information. In the present study, we want to use the molecular dynamics technique to obtain these parameters via computer simulations. Today’s advanced computer technology allows us to model simplified micelles, and a fair number of such studies can be found in the literature (e.g., refs 7-11). Nevertheless, atomistic molecular modeling on workstation computers is still limited to a few thousand atoms and to the nanosecond time scale. This usually prohibits the full simulation of micellar assembly for realistic molecules and, of course, the modeling of the phase behavior of such systems. An alternative possibility is to model a small portion of a micelle only using the proper boundary conditions to extract local information, which then can serve as input to, for instance, analytical theories. This is what we do here. The specific system which we study is 2,3,6,7,10,11-hexakis(1,4,7-trioxaoctyl)triphenylene (TP6) in water (cf. Figure 1). The monomer is a disklike molecule consisting of an aromatic core surrounded by hydrophilic side chains at its periphery. Above the critical micelle concentration, the monomers stack (the reason for this is one aspect of this work) and thus reversibly aggregate to form one-dimensional or rodlike micelles (see Figure 1). This system has been studied extenX
Abstract published in AdVance ACS Abstracts, July 1, 1996.
S0022-3654(95)03790-7 CCC: $12.00
Figure 1. 2,3,6,7,10,11-Hexakis(1,4,7-trioxaoctyl)triphenylene (TP6) molecule. Left: The molecular structure is drawn including the atom types corresponding to the notation used in Table 1. Upper right: Van der Waals representation of the TP6 molecule. The gray atoms are carbons, the dotted atoms denote oxygens, and the white atoms are hydrogens. Lower right: Sketch of the self-assembly process of TP6 molecules in aqueous solution. The stacking maximizes the contact between the cores of the TP6 molecules.
sively by Boden and co-workers12-17 (and references therein), which makes it an interesting model system for the above concept for combining theory and simulation. In the present work, we simulate a section of a onedimensional micelle representing a segment within an isolated, long micelle as well as a complete short micelle, both solvated by molecular water. We study the side-chain-water interface and deduce the effective diameter of the micelle. We also investigate the structure and the dynamics of the monomers within the micelle. In our simulation, the length of the simulated micellar segment is allowed to fluctuate freely subject to the external hydrostatic pressure. This allows us to determine the equilibrium monomer-monomer separation along the micelle’s axis, which is important for relating the aggregation number of a micelle to its length. Using a previously developed approach,18 we use the conformation statistics of the simulated segments to construct long micelles, for which we obtain the psersistence length for a number of different temperatures. Finally, we study the micelle size distribution via the free © 1996 American Chemical Society
MD Simulation of TP6 in Water
J. Phys. Chem., Vol. 100, No. 30, 1996 12163
enthalpy gain per monomer-monomer contact within a micelle using a thermodynamic integration technique. II. Methodology In molecular dynamics (MD) simulations, the Newtonian equations of motion
d2 1 b ri ) - ∇ B V (b r 1,...,b r N) (i ) 1, ..., N) mi dt2
(1)
are integrated numerically for all N atoms of the system located at positions b ri and having masses mi. Here we employ the program package AMBER 4.0 (Assisted Model Building with Energy Refinement)19 to compute the corresponding trajectory. The potentional V used in AMBER 4.0 consists of five different contributions; i.e.,
r N) ) V(b r 1,...,b
∑
(i) (i) 2 (i) 2 k(i) ∑ b (bi - b0 ) + ∑ kR (Ri - R0 ) + bonds i angles i
(i) (i) k(i) d [1 + cos(n φi - γ )] +
dihedrals i
∑
( )
+
∑
(2)
Aij
atom pairs ij
-
rij12
atom pairs ij
Bij
rij6 qiqj rij
The first two terms model all bond length bi and valence angle Ri deformations in terms of harmonic potentials. Here we use the SHAKE algorithm to constrain the bond lengths to their 20 The third term together with the 1-4 equilibrium values b(i) 0 . nonbonded interactions (cf. below) approximates the torsional potential variations in terms of the torsion angle φi. The fourth and fifth terms describe nonbonded interactions, i.e., LennardJones and Coulomb pair interactions. The summation over atom pairs includes all atom pairs separated by 3 (1-4 interactions) or more bonds in the same molecule or atom pairs where each bi - b rj|. atom belongs to a different molecule. Note that rij ) |r The nonbonded interactions are only calculated within a residuebased cutoff radius Rcut. This means that if two atoms belonging to two different residues are closer than Rcut, then all pair interactions between the two residues are calculated. In the following, each water and each TP6 molecule constitute one residue. Note also that, in addition, the 1-4 interactions are scaled by a factor of 1/2.19 The Lennard-Jones parameters Aij ) ijσij12 and Bij ) 2ijσij6 for the mixed interactions are obtained via the Lorentz-Berthelot mixing rules σij ) σi + σj (where a factor of 1/2 is absorbed into σi and σj) and ij ) (ij)1/2.21 The TP6 molecule is simulated using the united atom representation; i.e., CH, CH2, and CH3 groups are modeled by single pseudoatoms. Thus, 1 TP6 molecule consists of 66 atoms (i) (i) and pseudoatoms. The force-field parameters (k(i) b , b0 , k R , (i) (i) (i) (i) R0 , kφ , n , γ , σi, i) which we use for the TP6 molecule are taken from the OPLS/AMBER (optimized potentials for liquid simulations) united atom database22,23 as far as they are provided. It is worth noting that the original OPLS parametrization22 uses different mixing rules. However, in the present case, the resulting differences are less than 1% and, thus, can be neglected. The missing parameters are calculated using the semiempirical quantum chemistry package AM1.24 To do this, we first energy-minimize the molecular structure of a phenylmethyl ester molecule without constraints. This molecule is much simpler than TP6, but it contains all essential groups necessary for the determination of the missing parameters. Then we calculate the variation of the heat of formation as a function of the parameter of interest, keeping the rest of the molecule
TABLE 1: Force-Field Parameters Used in This Simulationa atom type
mass, amu
σ, Å
, kcal/mol
ref
C2 C3 CA CD OS HA OA
14.030 15.030 12.010 13.020 16.000 1.008 16.000
2.1916 2.1944 2.1046 2.1046 1.6837 0.0000 1.7766
0.1180 0.1600 0.1100 0.1100 0.1700 0.0000 0.1554
OPLS OPLS OPLS OPLS OPLS SPC/E SPC/E
bond
b30, Å
ref
CA-CA CA-CD CA-OS OS-C2 C2-C2 OS-C3 OA-HA HA-HA
1.400 1.400 1.380 1.425 1.526 1.425 1.000 1.633
OPLS OPLS AM1 OPLS OPLS OPLS SPC/E SPC/E
angle
kR, kcal/(mol rad2)
R0, deg
ref
CA-CA-CD CA-CD-CA CD-CA-OS CA-OS-C2 OS-C2-C2 C2-OS-C2 C2-OS-C3 HA-OA-HA HA-HA-OA
85.0 85.0 80.0 100.0 80.0 100.0 100.0
120.0 120.0 120.0 120.0 109.5 111.8 111.8 109.4 35.3
OPLS OPLS AM1 AM1 OPLS OPLS OPLS SPC/E SPC/E
dihedral
kn, kcal/mol
n
γ, deg
ref
X-CA-CA-X X-CD-CA-X X-CA-OS-X X-OS-C2-X X-C2-C2-X
1.325 2.650 8.000 1.450 2.000
4 2 1 1 1
180 180 90 0 0
OPLS OPLS AM1 OPLS OPLS
a Most of the parameters used for the TP6 molecule were taken from the OPLS/AMBER database.22,23 Parameters missing in the database are fitted via the method described in the text using the semiempirical quantum chemistry package AM1.24 The SPC/E water parameters are taken from ref 25.
fixed. From the resulting energy vs parameter curve, the missing parameter is obtained by fitting to the entire OPLS/ AMBER potential. The numerical values of all parameters are compiled in Table 1, and the corresponding atom types are shown in Figure 1. The force-field parameters used for the water molecules are taken from ref 25, and they are described in the next section. The partial charges of the atoms of the TP6 molecule are calculated with the charge equilibration algorithm provided with the POLYGRAF software package.26 Initially, the charges are obtained for the all-atom molecule, and subsequently, the charges of the hydrogens are transferred to the centers of the respective pseudoatom. In this fashion, we compute the charges for various conformationally different, separately minimized molecular structures and then average the results. This procedure together with the condition of overall neutrality yields a partial charge of qO ) -0.6e for each oxygen atom in the side chains and a charge of qC ) +0.3e for all carbon atoms bound to one of the side-chain oxygens. All other atoms of the TP6 molecules are uncharged. The partial charges of the water molecules are those of the SPC/E model discussed in the next section. We integrate the equations of motion (1) using the leap-frog verlet algorithm,21 and we apply periodic boundaries using the minimum image convention to calculate the nonbonded interac-
12164 J. Phys. Chem., Vol. 100, No. 30, 1996
Bast and Hentschke
tions. During the simulations, the temperature in the simulation box is controlled by the weak coupling velocity scaling thermostat of Berendsen et al.27 using a temperature relaxation time of 0.1 ps. In addition, the pressure is controlled by position scaling27 using a pressure relaxation time of 0.1 ps together with the default value of the inverse compressibility of 44.6 × 10-6 bar-1. First, the pressure is computed for each dimension separately, which is necessary due to the uniaxial anisotropy induced by the presence of the micelle. Then the atom positions are scaled for each dimension separately according to the respective pressures. This anisotropic pressure scaling is necessary because of two different effects. First, water molecules penetrate into the side chains of the micelle and the simulation box thus shrinks. With the pressure scaling perpendicular to the micellar axis, we achieve the proper water bulk density away from the micelle (as discussed below in the context of Figure 4). Second, the monomer-monomer separation is a priori unknown. By applying a constant pressure along the micelle, any stretch or compression of the micelle is allowed to relax, and the average monomer-monomer separation can fluctuate around its equilibrium value. III. Test of the Water Model We use the SPC/E water model developed by Berendsen et al.25 In this model, the OH bond length is 1 Å and the HOH valence angle is 109.47°. The hydrogen atom carries a partial charge qH ) 0.4238e. Correspondingly, the oxygen atom has a charge qO ) -0.8476e. The nonbonding Lennard-Jones interactions include interactions between the oxygen atoms only. The Lennard-Jones parameters are σOO ) 2 × 1.7766 Å and OO ) 0.1554 kcal/mol. Originally, the SPC/E water model was developed for Monte Carlo simulations assuming a rigid water molecule. Here, we introduce an additional bond between the two hydrogen atoms into the SPC/E model. By constraining this bond together with the two OH bonds using the SHAKE algorithm,20 the water molecule becomes rigid in our MD simulations too. This rigidity allows a longer time step during the simulation. The underlying assumption, of course, is that the behavior of the simulated system does not depend on the OH bond oscillations. The main advantages of the SPC/E model are its simplicity, which allows fast computation of the interactions, and its very good reproduction of experimental data for liquid water like the temperature dependence of the density and that of the self-diffusion coefficient as well as the O-O pair correlation function. Figure 2 summarizes the results of several 0.5-ns constantpressure test simulations of pure SPC/E water for different box sizes (i.e., 125 and 216 water molecules) at various temperatures ranging from 250 to 360 K. Figure 2a shows the simulated density as a function of temperature in comparison to the experiment.28 The deviation of the liquid densities from their experimental values is always less than 2%. For 216 water molecules and at T > 300 K, the deviation is even less than 1%. Also, for 216 water molecules, the figure includes the results for two different values of the cutoff radius, i.e. Rcut ) 7.5 and 9 Å. The reduction of the cutoff radius leads to a very minor decrease in density due to the neglect of attractive interactions. Figure 2b shows the corresponding temperature dependence of the self-diffusion coefficient, D, in comparison to the experiment.29 For liquid water, the self-diffusion coefficients are in very good agreement with the experimental data, especially in the temperature range which is of interest in this work (T ) 280-300 K). Notice that we do not observe any indication of freezing in the simulated temperature range below 273 K. Various reasons, like kinetic or finite size effects, are
Figure 2. Test of the SPC/E water model. Upper panel: Simulated water density vs temperature for different systems, i.e., 2 different volumes containing 125 (triangles) and 216 (circles and squares) water molecules, respectively. For the case of 125 molecules, the cutoff radius is Rcut ) 7.5 Å, whereas for the 216 molecules system, we compare two different cutoff radii; i.e., Rcut ) 9 Å (squares) and Rcut ) 7.5 Å (circles). The experimental values (diamonds) are taken from ref 28. Lower panel: Self-diffusion coefficients D of the water molecules vs temperature for two different box sizes, i.e., 125 molecules using Rcut ) 7.5 Å (triangles) and 216 molecules using Rcut ) 9 Å (squares). The experimental values (diamonds) are taken from ref 29.
thinkable, but their respective relevance is difficult to assess from the above simulations. However, it is worth mentioning that there have been simulations carried out in which the water molecules initially were arranged into an ice crystal structure. This structure remained stable during the simulation.30 The O-O pair correlation function, which is not shown here, is very similar to that obtained in Monte Carlo studies.31 The positions of the first three peaks of the experimental curve are reproduced to within 0.1 Å. Only the height of the first peak is somewhat exaggerated (roughly by 30%). In addition to the SPC/E model, we also investigated the TIP4P model of Jorgensen et al.31 A disadvantage of this model for MD simulations is that the additional charge site has to be treated as a fourth, massless atom. This leads to difficulties when integrating the equations of motion (1). In the literature, numerous other water models are discussed which include more detailed descriptions of the water molecule. For example, there have been attempts to simulate the lone pairs on the oxygen atom as additional charge sites,32 to introduce additional equations of motion for the partial charges,33 and to model the polarizability of the OH bond via additional dipoles along the bonds.34 However, in the present context, all these models are too complex and, thus, computationally too expensive. IV. The TP6/Water System For the various simulations in this work which include TP6, the starting configurations are prepared as follows. Stacks consisting of six and eight TP6 molecules are assembled by first building and energy-minimizing a single TP6 molecule with the program INSIGHTII using the Discover version 3.2 force field.35 The TP6 molecule is then replicated employing a 4-Å repeat distance perpendicular to the triphenylene plane. Subsequently, the stacks are surrounded by water molecules placed on a simple cubic lattice. The box with the 6-mer contains 2500
MD Simulation of TP6 in Water water molecules, whereas the box with the 8-mer contains 1827 water molecules. Note that both systems are prepared large enough to include bulk behavior of the water at large radial separations from the micelle. In the case of the 6-micelle, this is also true along the axial direction. Notice also that we apply a constant pressure of 1 bar to the box. In the case of the 8-micelle, the axial box dimension coincides with 8 times the monomer-monomer separation. The latter, however, is flexible due to the constant-pressure condition. Note that because of periodic boundary conditions, the 8-mer corresponds to a segment of a virtually infinite micelle. In a similar fashion, a third system is prepared which contains 2 isolated TP6 molecules solvated in 2004 water molecules. The results reported below are for MD simulations at three different temperatures (T ) 280, 290, and 300 K). After a short 10-ps MD run using a 1-fs time step, the remainder of each simulation is carried out with a 4-fs time step. The total length of each run is between 0.6 and 1.5 ns, during which the atom positions are stored every 0.2 ps. As an example, Figure 3 shows an instantaneous configuration during the simulation of the 8-mer. The eight monomers in the box form a stack with the core regions in contact. Notice that the water molecules penetrate into the side chains but do not enter the core region. Theoretical models2-6 for the phase behavior of rodlike micelles depend on model parameters like the micellar diameter, the monomer-monomer separation within a micelle, which determines its length, the micellar flexibility, and the effective monomer-monomer interaction within a micelle. Even though it is not yet possible to simulate the full phase diagram of an amphiphilic system like TP6 in water, the following simulation results can be used to obtain the above parameters. A. Side-Chain-Water Interface. In Figure 4, the radial density distributions of the water oxygen atoms and the radial distribution of the side-chain oxygen atoms are shown. Note that the water oxygen density is normalized to its experimental bulk value. The results plotted in Figure 4 are taken from the trajectory at T ) 300 K. The radial distributions extracted from the other two trajectories do not show significant deviations from the curves shown in Figure 4. The radial distribution function of the water molecules shows two peaks at distances of r ) 6 and 9.5 Å. In this range, the average density is approximately half that of the bulk value, which is approached in the range 11 Å < r < 14 Å. From this, we estimate a micelle diameter of about 28 Å. Notice that an appreciable onset of the water distribution occurs at r ) 4 Å; i.e., these water molecules are right at the periphery of the triphenylene core. However, the water molecules do not penetrate into the core region between the monomers. The corresponding radial distribution of the side-chain oxygen atoms shows three peaks which can be identified with the three oxygen atoms in each side chain. The first peak is sharp, whereas the broadness and the inward shift of the two outer peaks reflects the conformational freedom of the side chains. Notice that the arrows in Figure 4 denote the radial oxygen positions in a fully extended side chain. Comparing the two radial distributions, it can be seen that the radial density of the water oxygens shows local minima at the preferred positions of the side-chain oxygen distribution peaks. The oxygen atoms in the water molecule as well as in the side chains carry negative partial charges. Therefore, a repulsive Coulomb force acts between these atoms. This force explains the minima mentioned above. Nevertheless, there appear to exist hydrogen bonds between the neighboring water molecules and the side-chain oxygen atoms. This can be seen from the pronounced first peak at ≈1.9 Å in the side-chain
J. Phys. Chem., Vol. 100, No. 30, 1996 12165
Figure 3. Snapshot of the simulation box at T ) 300 K containing a stack of eight TP6 molecules in van der Waals representation surrounded by water molecules in stick representation. The upper panel shows a view along the micelle’s axis, whereas the lower panel shows a side view of the simulation box. Note that due to the application of periodic boundary conditions, we model a segment of eight monomers within a virtually infinite micelle.
oxygen-water hydrogen pair distribution function shown in the inset of Figure 4. B. Monomer-Monomer Stacking Distance. The monomer-monomer distance in the micellar aggregates is important in converting aggregation numbers into the actual lengths of the micelles. These lengths play an important role in the analytical description of the phase behavior of TP6 in aqueous solution. Here, we study two differently defined monomermonomer separations. In the first case, we determine the separation of the centers of mass dC of adjacent triphenylene cores. In the second case, we compute the projection of the previously determined separation dC onto the backbone in order to obtain the stacking distance dS along the backbone. In this study, we define the backbone of the micelle as the principal axis of inertia corresponding to the direction along the micelle. Figure 5 shows dC and dS averaged over the monomermonomer contacts within the periodic 8-micelle as a function of simulation time t for two temperatures, i.e., T ) 280 and
12166 J. Phys. Chem., Vol. 100, No. 30, 1996
Bast and Hentschke
Figure 4. Radial density distribution of the water oxygen atoms F (normalized to its experimental bulk value F0, solid line) and the radial distribution of the side-chain oxygen atoms ∆n/∆r (in arbitrary units, dotted line) at T ) 300 K. Note that r is the radial distance from the micelle’s backbone. The arrows indicate the radial positions of the oxygen atoms in a fully extended side chain. From the F/F0 curve, we estimate a radius of the micelle of 14 Å. The respective area of the plot is shaded. Inset: Side-chain oxygen-water hydrogen pair distribution function g(r′ ) (solid line) and the corresponding side-chain oxygen-water oxygen g(r′ ) (dotted line).
Figure 5. Upper panel: Center-of-mass separation dC between adjacent monomers in the 8-micelle vs simulation time t. Lower panel: Corresponding stacking distance dS vs simulation time t. The stacking distance is the center-of-mass separation projected onto the backbone of the micelle. Circles, T ) 280 K, squares, T ) 300 K. All values are obtained by averaging over all monomer-monomer contact sites within the micelle as well as over 100-ps time slices along the trajectory. The error bars denote the standard error assuming that the contact sites are independent of each other.
TABLE 2: Mean Values for the Center-of-Mass Separation dC and the Stacking Distance dS between Two Adjacent Monomers for the Three Different Temperatures T, K
〈dC〉, Å
〈dS〉, Å
280 290 300
4.5 ( 0.14 4.5 ( 0.19 4.4 ( 0.10
4.2 ( 0.12 4.3 ( 0.20 4.2 ( 0.13
300 K, which bracket the temperature range investigated here. The average over all contacts for the separations dC and dS are given in Table 2. On average, dC exceeds dS by 0.2-0.3 Å. However, we do not see any change of the mean distances with temperature exceeding the error margin. Thus, the magnitude of the error, where we assume that different contacts are statistically independent, is certainly an upper limit for the
Figure 6. Upper panel: Average center-of-mass separation dC between adjacent monomers for all eight individual contact sites n within the 8-micelle. Lower panel: Respective average stacking distance dS. The stacking distance is the projection of the center-of-mass separation onto the backbone of the micelle. Circles, T ) 280 K; hollow squares, T ) 290 K; solid squares, T ) 300 K. The error bars (only shown for T ) 290 K) denote 1 standard deviation.
temperature dependence of dC and dS. It is worth noting that the averages for dC and dS in Table 2 are rather close to the experimentally determined ring-ring separations of between ≈4.2 Å at T ) 280 K and 5 Å at T ) 300 K.17 In Figure 6, the two time-averaged separations for all eight individual contact sites within the periodic 8-micelle are shown for three temperatures T ) 280, 290, and 300 K. The errors (only shown for T ) 290 K) here denote 1 standard deviation. It is worth mentioning that the monomer-monomer separations appear to develop a superstructure not present in the starting configuration of the simulation; i.e., the monomer-monomer separation prefers to alternate between a “short” and a “long” distance. In addition, we compare the results obtained for the finite micelle with those for the micelle consisting of six monomers only. Within the fluctuations, we do not observe any difference in either dC and dS between the short and the finite micelle. This is also true if one compares the monomer-monomer separations at the two ends and in the middle of the 6-micelle. Thus, we conclude that in this case, our assumption that there is only a negligible difference between micelles of finite and infinite length holds. C. Intramicellar Rotational Diffusion. We are also interested in the dynamics and the mutual orientation of the monomers within the micelle. Here, we study the rotation diffusion of the monomers with respect to the long axis of the micelle and the distribution of the rotation angle between adjacent monomers. In order to study the rotation of the monomers within a micelle, we compute a one-dimensional rotation diffusion coefficient R, which we define in analogy to the translation selfdiffusion coefficient; i.e.,
R ) lim tf∞
〈(β(t) - β(0))2〉 2t
(3)
Here, the angle β(t) is defined in terms of a molecule-fixed axis in the plane of the monomer and a space-fixed axis (cf. sketch
MD Simulation of TP6 in Water
Figure 7. Mean squared rotation angle 〈∆β2(t)〉 of the monomers with respect to the long axis of the micelle as a function of time t. Thick solid line, T ) 280 K; dashed line, T ) 290 K; dotted line, T ) 300 K. The straight solid lines are linear fits to the simulation results for large t. Lower right: The sketch illustrates the definition of β(t) where the axis of rotation is defined by the principal axis of inertia of the 8-micelle corresponding to the micelle’s backbone. Upper left: Histogram of the angles δ between neighboring monomers (cf. the sketch). Dark shading, T ) 280 K; lighter shading, T ) 290 K; lightest shading, T ) 300 K.
in Figure 7). To only include the rotation around the micelle’s long axis, β(t) is computed between the projections of the two aforementioned axes onto a plane perpendicular to the micelle’s backbone. The angle β(t) thus is the net angle of rotation around the micelle’s backbone for a single monomer after time t. Note that there is no significant net rotation of the micelle as a whole with respect to the surrounding solvent. In Figure 7, we plot the mean squared rotation 〈∆β2〉, where ∆β ) β(t) - β(0), as a function of time t for T ) 280, 290, and 300 K. From the slope of the curves for large t, we determine R(T). For T ) 280 K, we obtain R ) 0.157 deg2/ps, for T ) 290 K, R ) 0.224 deg2/ps, and for T ) 300 K, R ) 0.318 deg2/ ps. As expected, with increasing T, the monomers start to rotate faster. Note again that although we have measured β with respect to a space-fixed frame, these values of R refer to the rotation of the monomers within the micelle, because the micelle itself does not rotate. Figure 7 also includes a histogram showing the distribution of the rotation angle δ of the monomers with respect to their neighbors. Here δ ) 0 means that the molecule-fixed axes of the two neighboring monomer have parallel projections onto a plane perpendicular to the micelle’s backbone (cf. sketch in Figure 7). It can be seen that the monomers prefer an eclipsed orientation which maximizes their core contact area rather than a staggered one. There appear to be differences between the three temperatures, but the limited statistics prohibit their interpretation. D. Micellar Flexibility. To obtain the flexibility of the rodlike micelle, we estimate its persistence length, P. P is defined by the correlation between the tangent vectors along the micelle (see upper inset of Figure 8). For not too small ui+n separations, the correlation for two tangent vectors b ui and b
J. Phys. Chem., Vol. 100, No. 30, 1996 12167
Figure 8. Orientation correlation function 〈u bi‚u bi+n〉i vs distance n. The b ui and b ui+n are tangent vectors along the micelle’s contour at monomers i and i + n. The average is taken over i. The smooth straight line is an exponential fit to the simulation result for T ) 300 K. Upper right insert: Definition of a tangent vector along the micelle’s contour. Lower left insert: x,y projection of the contour of one large micelle.
separated by n monomers is given by
〈u bi‚u bi+n〉i ) e-n/P
(4)
(cf. ref 36). In order to compute the persistence length via this equation, a large micelle is required to compute the average in (4). Here, we use an adaptation of a method which we have developed previously to compute P for a helical polypeptide (poly(γ-benzyl L-glutamate)).18 The method is a build-up procedure by which a long micelle is constructed by chaining together instantaneous conformations taken from a simulation of a short segment, i.e., the 8-micelle in the present case. In the first step, we extract the instantaneous configuration and attending conformations of stacks of four consecutive monomers of the 8-micelle. Note that we restrict the number to four out of a maximum of eight in order to avoid correlation effects due to the periodic boundary conditions. In the second step, we assemble these 4 segments into 1 large micelle consisting of 2000 monomers in the present case. The build-up procedure is based on finding the smoothest transitions from one segment to the next, i.e., the smoothest transition between the four segments to be added and the corresponding end segment of the micelle. As the new segment, we choose that segment from the list of previously stored segments whose first two monomers have the most similar triphenylene configuration compared to the upmost two monomers of the already existing micelle fragment. In addition, we impose one further condition. The temporal separation of two segments along the trajectory which during the buildup become consecutive in the micelle has to be at least 100 ps to avoid time correlations. Once the new segment is found, the two overlapping monomers are removed before adding the segment to the existing micelle fragment. Starting from a randomly chosen segment, the long micelle is thus assembled following the above procedure. The persis-
12168 J. Phys. Chem., Vol. 100, No. 30, 1996
Bast and Hentschke confinement must be included). The last term constitutes a simple one-parameter model of monomer-monomer interaction within a one-dimensional n-micelle, where n - 1 is the number of contacts between monomers. The quantity -R is the molar free enthalpy per monomer-monomer contact and RT, i.e.,
R)-
Figure 9. Average persistence length P for the constructed long micelles vs temperature T for three different temperatures.
∆Gcontact RT
which in general depends on the type of solvent, temperature, pressure, etc., and usually also on n, because monomers near the end of a micelle feel a different environment compared to those in the bulk of the micelle. Here, however, we will be exclusively dealing with the bulk value of R, which is independent of n. Using eq 5 in conjunction with the equilibrium condition µn ) nµ1 yields the micellar size distribution 0
tence length is then calculated using eq 4 by averaging the dot ui+n along the micelle. products of the tangent vectors b ui and b The direction of each tangent vector is defined as the moment of inertia axis along a cylindrical segment of the micelle consisting of five monomers centered at the origin of the tangent vector (cf. inset in the upper right of Figure 8). As a measure of the uncertaintiy of the persistence length to be computed, we calculate the standard deviation of P based on 10 different long micelles each constructed as described by choosing at random 10 different initial segments. Figure 8 shows an bi+n〉i, for 1 of the 10 micelles at T ) 300 exponential fit to 〈u bi‚u K. The lower left inset of Figure 8 shows an example micelle in terms of its contour’s x, y projection. In Figure 9, the mean value of the persistence length is shown for three different temperatures. It can be seen that the persistence length decreases with increasing temperature. Thus, the micelles become more flexible with rising temperature, which is consistent with the increase of the monomer rotational mobility within the micelle. The error bars denote the standard error as explained above. In this context, it is also worth noting that there is an additional effect due to the fact that the long micelle is assembled from 4-mers as explained above. The direct coupling between TP6 molecules separated by more than four monomers is therefore neglected. However, in comparison, the perhaps most widely used method for analyzing polymer conformations, the rotational isomeric state model,37 is usually based on the interactions between nearest neighbors only. E. Thermodynamics of Micellar Assembly. To good approximation, the chemical potential µn of a one-dimensional n-micelle in dilute solution can be written as
1 X + RT ln γ - RRT(n - 1) n n
( )
µn ) nµ˜ 0n + RT ln
(6)
0
Xn ) nXn1eR(n-1)+n(µ˜ 1-µ˜ n)/RT
(7)
under the given conditions. In the following, we wil be concerned with the calculation of R or rather ∆Gcontact, which, as eq 7 shows, is an important ingredient for the calculation of the micellar size distribution. Also, the analysis of the experimental data in ref 17 indicates that for the present system, the approximation µ˜ 01 ) µ˜ 0n seems to be justified and that the size distribution is mainly determined by R alone. In the following, we apply the simulated thermodynamic integration approach (cf. e.g., ref 38), which allows us to evaluate ∆Gcontact via the relation
∆Gcontact )
〈 〉
λeq ∂HC 1 dλ ∫ ∞ n-1 ∂λ λ
(8)
Here HC(λ) is the Hamiltonian of the system with a fixed monomer-monomer separation λ. 〈...〉λ denotes an ensemble average at the same fixed λ. The numerical integration of the right-hand side of eq 8 requires a series of MD simulations to evaluate the ensemble average at several distinct values of λ. However, the method is computationally too expensive for the complete system containing the micelle and the solvent. Therefore, we calculate ∆Gcontact along the alternative route (double lined arrows) in the following thermodynamic cycle n-micelle in solution with λ=∞
(n–1)∆Gcontact
hyd –n∆Gmon
n-micelle in solution with λ = λeq ∆Gnhyd –mic
(9)
(5)
(an extensive discussion of this equation can be found in the article by A. Ben-Shaul and W. M. Gelbart in ref 1). In the first term, µ˜ 0n is a standard average chemical potential per monomer in an n-micelle, excluding the interaction between the monomers, which is treated separately. The second term is the contribution due to the solute mixing entropy. Note that here Xn refers to the mole fraction of monomers within micelles of aggregation number n. In particular, n ) 1 corresponds to the free monomers. The third term describes the interaction between the solute particles. However, here we are merely interested in the limit of no interaction, where the activity coefficient γ ) 1 so that this term vanishes (at high concentration, however, excluded volume effects as well as entropic repulsion between the persistent-flexible micelles due to their
n-micelle in vacuum with λ=∞
vac (n–1)∆Gcontact
n-micelle in vacuum with λ = λeq
i.e., hyd vac (n - 1)∆Gcontact ) -n∆Gmon + (n - 1) ∆Gcontact + hyd (10) ∆Gn-mic hyd Here, ∆Ghyd mon and ∆Gn-mic are the free enthalpies of hydration of the monomer and the n-micelle, respectively, whereas vac ∆Gvac contact is ∆Gcontact in vacuum. Note that ∆Gcontact requires substantially less computational effort, because without solvent the number of interactions is greatly reduced.
MD Simulation of TP6 in Water
J. Phys. Chem., Vol. 100, No. 30, 1996 12169
Before we proceed, it is instructive to apply eq 10 under the following simplifying assumptions. First, we neglect the side chains completely. This assumption, which is discussed in detail below, implies that the (positive) difference - [n/(n hyd hyd + [1/(n - 1)]Gn-mic,sidechains is balanced by an 1)]Gmon,sidechains overall attraction between the side chains in the micelle (which yields a negative contribution to ∆Gcontact). In addition, we assume T ) 0 K when calculating ∆Gvac contact; i.e., we replace by the contact potential energy ∆Evac ∆Gvac contact contact. The latter is not a severe approximation, because after we have made the first assumption, we are left with only the rigid triphenylene core. Thus, eq 10 becomes
∑i∆G(i)e-∆G /RT - ∑i∆E(i)e-∆E /RT ∑i e-∆G /RT ∑i e-∆E /RT (i)
∆Ghyd )
(i)
(i)
(i)
(13)
Here, i distinguishes the possible conformations (of either the monomer or the micelle). ∆G(i) is the free enthalpy of just one fixed conformation i in solution, whereas ∆E(i) is the internal energy of that conformation in the gas phase. Scheraga and co-workers39,41 have proposed a model where ∆G(i) is written as
∆G(i) ) ∆E(i) + ∑∆ghyd k V(i)k
n hyd vac ∆Gcontact ≈ + ∆Econtact + ∆Gtriph n-1
(14)
k
1 hyd (11) ∆Gn-triph n-1 hyd is the free energy of hydration of the isolated Here, ∆Gtriph hyd is the free enthalpy of triphenylene molecule and ∆Gn-triph hydration of the “bare” triphenylene “micelle”. For large n, the last term can also be neglected, because the water accessible volume of the n-mer is small compared to the combined water accessible volume of the isolated monomers. Thus, we finally obtain
hyd vac + ∆Econtact ∆Gcontact ≈ -∆Gtriph
free enthalpy of hydration:
(12)
hyd , we estimate based on the experiThe first term, - ∆Gtriph hyd mental value for ∆Gpyrene ) -4.46 kcal/mol taken from Table III of ref 39 (cf. also Table VII in ref 40 on which Table III in ref 39 is based). In order to estimate ∆Evac contact, we assume that only the interactions between nearest and second-nearest triphenylene rings contribute. The nearest-neighbor potential energy we obtain by energy-minimizing (here we employ the MM2 force field combined with a π-orbital SCF calculation to adjust the force-field parameters (CSC Chem 3D, Version 3.1, standard parameter settings, Cambridge Scientific Computing, Inc.) for higher accuracy) an eclipsed triphenylene dimer, which yields -37.7 kcal/mol. From this, we must subtract the self-energy of the individual triphenylenes, which is -12.5 kcal/mol. Thus, the dimer contact energy is ∆Evac contact(dimer) ) -12.7 kcal/mol. It should be mentioned that the configuration of the energyminimized dimer closely resembles the stacking of adjacent planes in graphite. The energy contribution due to the nextnearest-neighbor interaction we estimate from an analogus minimization of a trimer, for which we obtain ∆Evac contact (trimer) ) -26.9 kcal/mol. The next-nearest-neighbor energy is then -26.9 - (2 × -12.7) kcal/mol ) -1.5 kcal/mol. With this, we finally have ∆Evac contact ≈ -12.7 + (-1.5) ≈ -14.2 kcal/mol and ∆Gcontact ≈ 4.5 - 14.2 ) -9.7 kcal/mol or RT)300K ) 16.3 (using RT ) 0.596 kcal/mol at T ) 300 K). This is remarkably close to the corresponding experimental value R ≈ 13.8.17 However, in the next sections, we find that the neglected contributions are rather large and do not completely cancel. Thus, the close accord of this simplified calculation with the experiment must be considered accidental. hyd hyd . In this section, we and ∆Gn-mic 1. Calculation of ∆Gmon use the hydration shell model introduced by Scheraga and coworkers39,41 together with our simulation trajectories to compute the hydration free enthalpies of an n-micelle and of a monomer. Their calculation is based on the following expression for the
i.e., the hydration contribution to ∆G(i) is written in terms of is a free enthalpy of atomic increments k, where ∆ghyd k hydration density per atom (or united atom) times its wateraccessible volume, V(i)k, in the conformation i. V(i)k is defined as the volume of the spherical hydration shell of thickness hyd - RvdW is its hydration shell radius and Rhyd k k , where Rk vdW Rk is the van der Waals radius of atom k, minus the overlap volume due to the other van der Waals spheres intersecting the hydration shell. The overlap and thus V(i)k depend on the conformation i. Scheraga and co-workers41,39 have tabulated and the attending values of Rhyd and RvdW for a large ∆ghyd k k k number of atomic groups by fitting the theoretical values of ∆Ghyd to the experimental ones using a set of training compounds. In the following, we use their values for ∆ghyd k , vdW as listed in Table I of ref 39. The conformaRhyd k , and Rk tion-dependent water-accessible volume V(i)k, however, we calculate with a Monte Carlo integration instead of their analytic approach. This is done by selecting a large number of points (here we use 100 000 points) in the simulation box at random. Then we determine the fraction of these points which belong to the hydration shell of atom i but do not simultaneously lie inside a van der Waals shell of any other atom of the micelle. By multiplying this fraction with the known volume of the simulation box, we obtain V(i)k. In order to obtain ∆Ghyd from (13), we make the simplyifying assumption that the weighting of the side-chain conformations in solution is the same as in vacuum. With this assumption, (13) and (14) yield
∆Ghyd ≈ ∑∆ghyd k 〈V(i)k〉
(15)
k
where 〈V(i)k〉 can be calculated directly based on our simulation trajectories above. Note again that ∆ghyd k is independent of the may monomer’s or micelle’s conformation. However, ∆ghyd k be altered due to the polarization of the hydration shell of k by goes over into ∆ghyd + a nearby polar group l. Thus, ∆ghyd k k hyd,pol hyd,pol is given by eqs 5-8 in ref 39. ∑l∆gk,l , where ∆gk,l depends on the separation rkl between k and Because ∆ghyd,pol k,l l and, therefore, on the conformation, (15) becomes
∆Ghyd ) ∆Ghyd,no-pol + ∆Ghyd,pol ) ∑ ∆ghyd k 〈V(i)k〉 + k
hyd,pol 〈∆gk,l V(i)k〉 ∑ k,l
(16)
To compute ∆Ghyd mon of the separated monomer, three additional simulations are carried out. During these runs, the simulation box contains only 2 TP6 molecules solvated in 2004 water
12170 J. Phys. Chem., Vol. 100, No. 30, 1996
Bast and Hentschke
vac Figure 10. Time-averaged contact force 〈fax 〉λ as a function of λ for three different temperatures T ) 280 K (squares), 290 K (diamonds), and 300 K (circles). For comparison, the result of a similar calculation for T ) 300 K only with a different charge model (AM1 charges) is also shown (crosses). All points are averages over MD trajectories of 0.5 ns.
TABLE 3: Free Enthalpies Calculated along the Alternate Path in (9) and ∆Gcontact and r Calculated According to (10) and (6), Respectively, in the Limit for Large n ∆G, kcal/mol hyd ∆Gmon
hyd (1/n)∆Gn-mic vac ∆Gcontact ∆Gcontact R
T ) 280 K
T ) 290 K
T ) 300 K
-47.8 -29.3 -43 -24.5 43.9
-47.6 -29.4 -38 -19.8 34.3
-48.4 -29.8 -35 -16.4 27.4
molecules. The separation between the monomers is always large enough so that their interaction is negligible. All other parameters are the same as in the original simulations. Before we discuss the results obtained via (16), however, we turn to the calculation of ∆Gvac contact. vac . Note that the right-hand side 2. Calculation of ∆Gcontact of eq 8 is the reversible work necessary to separate a monomer from the micelle. Here, we evaluate the right-hand side of (8) in vacuum by calculating vac ∆Gcontact ) ∫∞ 〈fvac ax 〉λ dλ λeq
(17)
where 〈fvac ax 〉λ is the time-averaged force along the axis of the micelle in vacuum acting on a monomer, which is pulled off the end of the micelle. The micelle is approximated by a trimer for which we perform molecular dynamics simulations (as discussed above) constraining the positions of the carbon atoms of the central hexagon in each monomer so that these hexagons in neighboring monomers are eclipsed and parallel. All other atoms of the triphenylene core and in particular the side chains remain unconstrained. Initially, we set the distance of the triphenylene to the equilibrium separations discussed above. Subsequently, the monomer is pulled away along the axis of the micelle, and 〈fvac ax 〉λ is calculated as a function of λ. hyd hyd vac , ∆Gn-mic , and ∆Gcontact . Table 3 3. Results for ∆Gmon lists the hydration contribution to ∆Gcontact according to eq 15 neglecting the polarization contribution (cf. eq 16). Thus, in the indicated temperature range, the hydration contribution to hyd ∆Gcontact, - ∆Ghyd mon + (1/n)∆Gn-mic, is close to 18.5 kcal/mol, which also means that the hydration disfavors micelle formation. If we compare this with our above approximate calculation (cf. eq 10), we can estimate that the side-chain contribution to this ≈14 kcal/mol. Table 3 also contains the values of ∆Gvac contact calculated according to (17) (the corresponding 〈fvac ax 〉λ vs λ curves are shown in Figure 10). The values for ∆Gvac contact obtained in this fashion range from -43 kcal/mol for T ) 280 K to -35 kcal/mol for T ) 300 K. Note that the decrease in
∆Gvac contact with decreasing temperature is physically reasonable, because the entropic repulsion (due to spatial confinement) between the side chains should become less as temperature decreases. If we take the -14.2 kcal/mol obtained above for ∆Evac contact as the core contribution, then the contribution due to direct side-chain interaction is between about -29 and -21 kcal/ mol. Thus, the magnitude of the side-chain hydration contribution to ∆Gcontact of about 14 kcal/mol is in fact considerably smaller than the magnitude of the direct side-chain interaction contribution. In particular, we see that it is the strong negative contribution of the direct monomer-monomer interaction which finally favors the micelle formation. The resulting values for ∆Gcontact yield an R (cf. Table 3) which is between 43.9 for T ) 280 K and 27.4 for T ) 300 K, and thus, our R is about 2-3 times larger than the experimental value quoted above. One source of error in the present calculation is the neglect of the polarization contribution to ∆Ghyd in eq 16. The parametrization is strongly model dependent and unfortunately cannot be carried over to the model used in this work. To estimate the corresponding error, we have calculated ∆Ghyd for some of the training compounds used for parametrization in ref 39. For small n-alkanes, our values coincide with the experimental values listed in ref 39 to within 1%. On the other hand, for 1,2-dimethoxyethane, which is similar to the monomer’s side chains, our value for ∆Ghyd exceeds the experimental value by 10%. Note, however, that in general, for polyfunctional molecules of N and O atoms, the discrepancy between the above method and experiment can be as large as about 40% (cf. Table XI of ref 39). We also tested the sensitivity of 〈fvac ax 〉λ to the charge model, because the partial charges are a possible source of considerable error in every force-field calculation. Figure 10 shows the integrand on the right-hand side of (17) for T ) 280, 290, and 300 K. All three curves show a similar behavior with a positive (repulsive) maximum at λ ≈ 7 Å. We have tested the sensitivity of 〈fvac ax 〉λ on the partial charges by comparing the result for the partial charges of the Q-equilibration algorithm, which we use here, to the result for AM1 charges (which on average are 2.6 times smaller). Overall, we find little difference in the 〈fvac ax 〉λ vs λ curves, and thus, the corresponding areas under the curves are close, and we would estimate the error at about 10%. Thus, hyd vac if the error for ∆Ghyd mon, (1/n)∆Gn-mic, and ∆Gcontact in Table 3 is 10%, then the resulting error for R is already about 25%, because ∆Gvac contact essentially is calculated as the difference between significantly larger numbers. V. Conclusion In this paper, we have studied the structure, dynamics, and thermodynamics of a rodlike micelle using molecular dynamics simulations. From the radial water density profile and the radial positions of the side-chain atoms, we estimate the effective diameter of the TP6 micelle in water to be 28 Å. For the average monomer-monomer center of mass and stacking distance within the micelle, we obtain 4.5 and 4.2 Å, respectively, for all three temperatures considered here. These average values are in good agreement with the experimentally determined ring-ring separations. In addition, we observe a superstructure characterized by alternatingly long and short average monomer separations along the micelle. Long TP6 micelles are flexible, and we obtain persistence lengths between ≈340 monomers at 280 K and ≈160 monomers at 300 K. It is worth mentioning that the temperature dependence of the phase behavior of micellar systems like the one studied here enters mainly through two quantitiessthe persistence length and the
MD Simulation of TP6 in Water contact free enthalpy. Here, we have shown that the temperature dependence of the first is rather strong in the case of TP6 in water; i.e., P is reduced by a factor of ≈2 between T ) 280 and 300 K. Thus, the decrease of the persistence length with increasing temperature is substantial. Also, the decrease of the persistence length is accompanied by an increased rotational mobility of the monomers within the micelle. As mentioned before, the effective micelle diameter as well as the persistence length are key quantities in the theoretical description of the phase behavior of the TP6/H2O system; i.e., excluded volume interactions are dependent on the micelle’s geometry (aspect ratio), whereas the orientation distribution of the micelle’s contour and the attendant entropy contribution depend on the persistence length and on the micelle’s contour length, which in turn depends on the monomer-monomer separation. The present simulation shows how, even for a large system (in terms of the monomer size), these quantities can be calculated including microscopic detail. Finally, we have considered the monomer-monomer contact free enthalpy which largely governs the micellar size distribution at low concentrations. We base this calculation on a combination of our simulation trajectories with a hydration shell model due to Scheraga and co-workers. Even though our simulation results of R ) 27.443.9 overestimate the experimental findings, they allow us to compare the various contributions to the monomer-monomer contact free enthalpy. We find that the main contribution to R is due to the competition between the monomer contact potential energy (which favors micellation) and the loss of free enthalpy of hydration when a monomer-monomer contact is formed (which disfavors micellation). References and Notes (1) Gelbart, W. M., Ben-Shaul, A., Roux, D., Eds. Micelles, Membranes, Microemulsions and Monolayers; Springer: New York, 1994. (2) Taylor, M. P.; Herzfeld, J. J. Phys.: Condens. Matter 1993, 5, 2651. (3) McMullen, W. E.; Gelbart, W. M.; Ben-Shaul, A. J. Chem. Phys. 1985, 82, 5616. (4) Odijk, T. J. Phys. 1987, 48, 125. (5) Hentschke, R. Liq. Cryst. 1991, 10, 691. (6) van der Schoot, P.; Cates, M. E. Europhys. Lett. 1994, 25, 515. (7) Bo¨cker, J.; Brickmann, J.; Bopp, P. J. Phys. Chem. 1994, 98, 712. (8) Rector, D. R.; van Swol, F.; Henderson, J. R. Mol. Phys. 1994, 82, 1009. (9) Laaksonen, L.; Rosenholm, J. B. Chem. Phys. Lett. 1993, 215, 429. (10) Shelley, J. C.; Sprik, M.; Klein, M. L. Langmuir 1992, 9, 916. (11) Karaborni, S.; Esselink, K.; Hilbers, P. A. J.; Smit, B.; Kartha¨user, J.; van Os, N. M.; Zana, R. Science 1994, 266, 254. (12) Boden, N.; Bushby, R. J.; Hardy, C. J. Phys. Lett. Paris 1985, 46, L 325.
J. Phys. Chem., Vol. 100, No. 30, 1996 12171 (13) Boden, N.; Bushby, R. J.; Hardy, C.; Sixl, F. Chem. Phys. Lett. 1986, 123, 359. (14) Boden, N.; Bushby, R. J.; Jolly, K. W.; Holmes, M.; Sixl, F. Mol. Crystallogr. Liq. Crystallogr. 1987, 152, 37. (15) Ferris, L. M. Ph.D. Thesis, University of Leeds, 1989. (16) Edwards, P. J. B. Ph.D. Thesis, University of Leeds, 1993. (17) Hubbard, J.; Boden, N. In preparation, 1995. (18) Helfrich, J.; Hentschke, R.; Apel, U. M. Macromolecules 1994, 27, 472. (19) Amber 4.0: Pearlman, D. A.; Case, D. A.; Caldwell, J. C.; Seibel, G. L.; Singh, U. C.; Weiner, P. K.; Kollman, P. A. Molecular Dynamics Simulation Package, University of California, San Francisco, 1991. (20) Ryckaert, J. P.; Cicotti, G.; Berendsen, H. J. C. J. Comput. Phys. 1977, 23, 327. (21) Allen, M. P.; Tildesley, D. J. Computer Simulations of Liquids; Oxford University: Oxford, 1987. (22) Jorgensen, W. L.; Tirado-Rives, J. J. Am. Chem. Soc. 1988, 110, 1657. (23) 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. (24) Dewar, J. D.; Zoebisch, E. G.; Healy, E. F.; Stewart, J. J. P. J. Am. Chem. Soc. 1985, 107, 3902. (25) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269. (26) Polygraf, Version 3.0. Molecular Dynamics Simulation Software, Molecular Simulations Inc., 1992. (27) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684. (28) Riddick, J. A.; Bunger, W. B.; Sakano, T. K. Organic SolVents; John Wiley & Sons: New York, 1986. (29) Scha¨fer, K., Ed. Landolt-Bo¨ rnstein, Vol. 5: Eigenschaften der Materie in ihren Aggregatzusta¨ nden; Springer: Berlin, Heidelberg, New York, 1969. (30) Karim, O. A.; Haymet, A. D. J. Chem. Phys. Lett. 1987, 138, 531. (31) Jorgensen, W.; Chandrasekhar, J.; Madura, J.; Impey, R.; Klein, M. J. Chem. Phys. 1983, 79, 926. (32) Stillinger, F. H.; Rahman, A. J. Chem. Phys. 1974, 60, 1545. (33) Rick, S. W.; Stuart, S. J.; Berne, B. J. J. Chem. Phys. 1994, 101, 6141. (34) Zhu, S.; Singh, S.; Robinson, G. W. J. Chem. Phys. 1991, 95, 2791. (35) Insight, Version 2.0, Molecular Dynamics Simulation Software, Biosym Technologies, 1989. (36) Landau, L. D.; Lifshitz, E. M. Statistische Physik, Teil 1; Akademie Verlag: Berlin, 1979; Vol. 5. (37) Mattice, W. L.; Suter, U. W., Eds. Conformational Theory of Large Molecules; John Wiley & Sons: New York, 1994. (38) van Gunsteren, W. F.; Berendsen, H. J. C. J. Comp.-Aid. Mol. Des. 1987, 1, 171. (39) Kang, Y. K.; Nemethy, G.; Scheraga, H. A. J. Phys. Chem. 1987, 91, 4110. (40) Cabani, S.; Gianni, P.; Mollica, V.; Lepori, L. J. Solution Chem. 1981, 10, 563. (41) Kang, Y. K.; Nemethy, G.; Scheraga, H. A. J. Phys. Chem. 1987, 91, 4105.
JP953790L
