J. Phys. Chem. B 2006, 110, 22719-22725
22719
Detailed Structures and Mechanism of Polymer Solvation Qing Ji and Xiaozhen Yang* Polymer Physics Laboratory, National Laboratory of Molecular Science, Institute of Chemistry, Chinese Academy of Sciences, Beijing 100080, China ReceiVed: July 9, 2006; In Final Form: September 6, 2006
In the present study, we simulated a model system, PE in biphenyl, to gain the insight into the detailed solvation structures and the molecular mechanism of polymer chain solvation. Using atomistic molecular dynamics (MD) simulation, it was found that when the biphenyl is far from PE chain or in the bulk, the dihedral angle of the two rings in the solvent molecule are ∼32°. But, the dihedral angel is about 27° when the biphenyls are very close to the PE chain. In the first solvation shell, the orientation angle of the biphenyl long axis to the chain segment backbone was found to be enhanced around two values: ∼0 and ∼60°. The detailed solvation structures found here include all dyad conformations (TT, TG, TG′, GT, GG, GG′, G′T, G′G, and G′G′) and vary as a function of the distance between PE chain and biphenyls in the first solvation shell. The closer the the solvent molecule to the PE segment, the higher the TT conformation fraction response is. The other dyad conformations such as TG, GG′, etc. undergo different decreases, respectively. This study shows that the solvation even in the Θ condition makes the overall size expansion or the chain stretched. Such a cooperative change was examined here and found not due to generating or losing a conformational state but due to a change in conformational distribution. This change occurs in the middle location of the chain instead of the chain end locations.
Introduction To our knowledge, solvation is a physicochemical phenomenon that a kind of molecule A (solute) evenly distributes into another kind of molecule B (solvent). The interaction between A and A finally occurs with some interactions between A and B. It is a process resulting in a solvation structure. For nonpolar systems, van der Waals forces predominate the intermolecular interaction, and the electronic structure in the ground state remains almost unchangeable. The process is a physical one. When a hydrogen bond or ion bond enters, the solvation structure should have slight or serious changes in its electronic structure. This change usually influences the chemical property of the solute molecule. The solvation becomes an impure physical process. Thermodynamically, the solvation results in a structural state, where the solutes interact with solvents. Statistically, there are thousands of interaction configurations in the solvation. Each of the configurations is called a solvation structure. Nowadays, people gradually perceive that the solvation structures, especially detailed solvation structures, are very important and are found in many disciplines, such as life science,1-4 pharmacy,5-8 lavation,9-13 etc. For example,4 the conformational stability of human polyclonal IgG in solution decreases with a decreasing pH; to know the reason, one must study the detailed solvation structure. The solvation structures are very helpful in designing medicine6,7 and detergents12,13 because the pharmacy molecule interacting with protein and detergent reacting with dirt are actually realized by special interaction configurations in the solvation. The polymer solvation structure has triggered numbers of studies for years.14 The pioneering works by Flory15 brought the thermodynamic theory of a polymer solution with a lattice * To whom all correspondences should be addressed. E-mail: yangx@ iccas.ac.cn.
model. This theory describes polymer chains dissolved in a solution. Generally, in the polymer solvation, people first focused their attention on a change in the overall size of the polymer chains.16-19 In this concept, the solvation structures refer to not only the overall size but also to the local conformation of both the polymer chain and the solvent, the interaction configuration between the solute and the solvent, and even the bond length and the bond angle. To discover detailed structures in polymer solution is not convenient in experiment and in theory. Experimentally, the overall size of the polymer chains can be measured by light scattering and neutron scattering. The local conformations can be, in principle, observed by using NMR,20,21 X-ray diffraction,22 and vibrational spectroscopy.23,24 Actually, one can obtain effective conformational data of polymer chains in solution by using NMR20,21 but cannot have signals with the same high resolution for the sample in the solid state. In the crystalline state, one can measure the conformational structures in the lattice by using X-rays. Using the same tool, one can also measure the conformational distribution of a molecule in liquid, combined with the spectral simulation by modeling the molecular conformations.22 As is well-known, FTIR23 and Raman24 spectroscopy are able to observe variations in polymer conformations since the overall size of the polymer chain is formed by the conformation of each σ-bond along the chain. Monitoring the conformations of all the backbone bonds along the chain is thus important. However, it is difficult to observe the conformation of every backbone bond in a repeat unit of a polymer chain. For example, poly(ethylene terephthalate) (PET) is one of the frequently used samples in conformational analysis of polymer chain by using infrared spectra.23 But, as there are seven σ-bonds in a repeat, one can only obtain conformational changes from one σ-bond, C-C, in ethylene glycol. Up to now, there was still no way to obtain quantitative analysis for the other
10.1021/jp0643101 CCC: $33.50 © 2006 American Chemical Society Published on Web 10/13/2006
22720 J. Phys. Chem. B, Vol. 110, No. 45, 2006 six σ-bonds. Fortunately, conformational variations of all the three σ-bonds in the poly(ethylene oxide) (PEO) repeat unit were observed24 by using the isotropic Raman spectroscopy and the method of the disordered longitudinal acoustic mode (D-LAM). This study first revealed that PEO solvation in water leads to a variation in conformational distribution of the three σ-bonds along the chains from the one in bulk. The difficulty in discovering the polymer solvation is in revealing the detailed structures, which explicitly describe the interaction configuration between polymer and solvent and have not been observed in experiments so far. In the theoretical works (perturbation theory,25 self-consistent field theory,26,27 and scaling theory, etc.28), the polymer solution was treated as a low density gas of polymers in a vacuum or in a continuum fluid. As Chandler and Pratt29 have pointed out, only those theories that explicitly consider the molecular nature of the solvent can hope to describe the interaction configurations; however, the fact is that at liquid densities, the environment of the solute molecule is highly structured, and without knowledge of the detailed structures, it is difficult to make a reasonable approximation. As a bridge linking experiment and theory, molecular simulations can reproduce most atomistic information about polymer solvation. It enables us to investigate the solvation structures in detail. Tasaki analyzed MD data30 and found three interaction configurations in which the mass center of the solvent molecule benzene is close to the PEO oxygen, and the distance is 4.1, 5.6, and 6.3 Å when the PEO segment is in different conformations, respectively. The discovery of such detailed structures was focused on in the present study as the first point. The second point we aimed at in the present study was the solvation mechanism. It is well-known that local conformations of polymer chains relate to the statistical overall size of the chains.17-19 When the polymer chains are in different states (i.e., the molten state, a solution with a poor solvent, or tahe solution with a good solvent), the local conformation changes. We think that such a conformation change is at least attributed to two mechanisms. One is to gain some new or lose some old conformational states (from t, g, g′ to t, t′, g, g′ or to t, g′), and the other is a change only in the conformational distribution without a change in the conformational state. To validate the two mechanisms, we must investigate the interaction structure between the polymer and the solvents in detail. In the present study, we built a dilute solution model of polyethylene (PE) in biphenyl and also used the model of molten PE as a comparison. The models underwent atomistic MD simulations. After validation of the force field for our system, detailed interaction configurations were analyzed by recognition of the first solvation shell of the solvent molecules to the polymer and also by the orientation distribution between the local segment and the solvent molecule. Furthermore, the features of the conformational sequences on various detailed interaction configurations were explored. Finally, we discussed the solvation mechanism of the PE chain in biphenyls. Models and Details of Simulation We built two cubic periodic boundary models, which stand for the PE solution and melts, respectively. The model for the solution was composed of one C45H92 molecule and 400 C12H10 (biphenyl) molecules. The model for the PE melts was composed of 45 C45H92 molecules. The cell length of the former is about 50 Å. That of the latter is about 42 Å. Both models were prepared for atomic MD simulations.
Ji and Yang The MD simulations were performed by using Material Studio version 3.2. The force field parameters and functions came from COMPASS31 (condensed-phase optimized molecular potential for atomistic simulation studies). Bonded potential terms involve the bond, the bond angle, and the torsion potential. Nonbonding interaction terms involve van der Waals and electrostatic interactions. The cutoff distance was set to 9 Å for van der Waals terms and 13 Å for electrostatic terms. The MD integration time step was 1 fs. The simulated temperatures were set at 400 K for the Θ solution of PE in biphenyl32 and 450 K for the molten state. These temperatures were controlled with the Nose´-Hoover method.33-35 The simulated pressures were at 1 atm. It was controlled with Andersen method.36,37 The recorded interval was set to 100 fs. To obtain equilibrium data, we ran a constant NPT simulation of 0.5 ns and a constant NVT simulation of 1.5 ns sequentially. In the subsequent simulation, the canonical simulation of 3 ns was run for sampling. Analyzed results are reported in the following sections. Results and Discussion Validation of the Results. Before beginning our discussion on the detailed structure and solvation mechanism of polyethylene in a biphenyl solution, we have to check the quality of the obtained data from the simulation. The checking includes the rationality of the nonbonding interaction and the duration length of the MD runs. Rationality of Nonbonding Interactions. The present simulation was performed at atomistic levels under classical mechanics, and therefore, using a reasonable force field is very fundamental. As of late, COMPASS31 has been considered as a rare but excellent force field. This is because COMPASS possesses carefully optimized nonbonding interaction parameters. This force field is different from general force fields, which pick up the nonbonding interaction parameters by either an arithmetical average or a geometrical one. In the present study, we are interested in intermolecular interactions between polymer solute and solvent, with regards to the very nonbonding interaction parameters. We therefore adopted COMPASS in the present work. As expected, the simulated results here appear reliable in many aspects, such as the density of the biphenyl liquid and the overall size of the PE chain. Biphenyl Liquid Density. Use of a reliable force field in building a molecular system results in a reasonable density, which is close to experimental observation. The density for room temperature was proven early.31 Since 400 K is our simulation temperature, it is also the Θ temperature of PE in the biphenyl solution.32 Therefore, it is necessary for us to examine the rationality of COMPASS in mimicking the biphenyl at 400 K. For this purpose, we built a bulk model of biphenyls. A total of 80 biphenyl molecules were placed in a cubic box with a three-dimensional periodic boundary condition. After a NPT run of 500 ps, the density of biphenyl was analyzed. The calculated density was 0.95 g/cm3 at 400 K. It is the same value as the experimental density of biphenyl bulk.32 It indicated that the COMPASS force field accurately depicts the nonbonding interaction between biphenyl molecules at 400 K. Overall Size of PE Chain. Usually, the overall size of a polymer chain is represented by the end-to-end distance or the radius of gyration. It is one of the basic properties that is directly influenced by nonbonding interactions. Here, we examined the overall size of the PE chain in biphenyl and in the molten state to examine the rationality of COMPASS in describing the nonbonding interactions between PE chains and between the polymer and the solvent molecule.
Structures and Mechanism of Polymer Solvation
J. Phys. Chem. B, Vol. 110, No. 45, 2006 22721
After analyzing the data from canonical simulations of 3 ns, we computed the statistical end-to-end distance, 〈R2〉0.5, and the statistical radius of gyration, 〈Rg2〉0.5. The 〈R2〉0.5 was calculated as 26.3 ( 3.5 and 27.2 ( 12.5 Å for the PE chain in the melt and in the biphenyl solution, respectively. The 〈Rg2〉0.5 was 10.2 ( 0.5 and 10.3 ( 3.5 Å. To validate the data, we examined 〈R2〉 based on the number of bonds n, the bond length l, and C∞ ) 6.8,38 where 〈R2〉 ) C∞nl2 ) 26.6. The value appears reasonable. The overall size value of PE in solution is bigger than that in melt in relation to the expansion factor, which measures the overall size of a polymer coil under special conditions reduced by the molten state. It is defined as
R2 ) 〈Rg〉s2/〈Rg〉m2
(1)
where the subscript signs s or m denote that Rg is measured in the solution or in the molten state. In general, when the polymer is in good solution, R > 1, when the polymer is in poor solution, R < 1, and when polymer is in an unperturbed solution, R ) 1. The experimental expansion factor of PE in biphenyl in the Θ condition ranged from 1.00 to ∼1.05.32 The present simulation was under the same Θ condition, and the obtained expansion factor is 1.01, which agrees with the experimental value. On the basis of the previous examinations, we concluded that COMPASS successfully describes the nonbonding interaction between PE chains and biphenyl molecules and that between PE and biphenyl. It is an appropriate force field for the present study. Effective Duration Length of MD Run. After confirming the force field, ascertainment turns to the duration length of the MD run. The run needs a sufficient period of time to have data equilibrated. Regularly, such a reasonable duration length can be estimated by using the time autocorrelation function (ACF).39 The ACF describes, during dynamics, how quickly a property portrayed by ACF forgets its initial value. A series of ACF were investigated here to ascertain whether the duration length of the MD run is sufficient in the present simulation. The components are the ACF for the dihedral angle φ, the radius of gyration Rg, and the end-to-end vector (r) as defined next.
Cφ(t) )
- 2 - 2
CRg(t) )
〈(Rg(t) - 〈Rg〉)(Rg(0) - 〈Rg〉)〉 (Rg(0) - 〈Rg〉)2 Cr(t) )
(2)
(3)
2
(4)
All decays of the three kinds of ACF for PE in the solution and in the melt are shown in Figure 1. This figure qualitatively shows that torsion φ needs the shortest time to reach the equilibrium, while r needs the longest time to relax. This result is in accordance with the results of Paul and Smith’s work.40 One can find that the PE chain experienced a longer relaxation time in the solution than in the melt, which indicates a damping effect from the smaller solvent molecules.30 To obtain quantitative relaxation time for φ, r, and Rg in the melt and in the solution, the Kohlrausch-Williams-Watts
Figure 1. Autocorrelation function of torsion, r, and Rg in the melt: (A) torsion, (B) Rg, and (C) r. In the solution: (D) torsion, (E) Rg, and (F) r.
TABLE 1: Averaged KWW Parameters and Relaxation Time for Torsion, r, and Rg of PE in Melt and Solution torsion A B τ (ps)
r
Rg
melt
solution
melt
solution
melt
solution
7 0.8 8
15 0.8 17
681 0.7 862
852 0.8 967
83 0.8 94
230 0.8 261
(KWW)41 function was used. The KWW function is defined as eqs 5 and 6
C(t) ) exp τ)
(-tR )
β
∫0∞ exp((- t/R)β)dt ) RβΓ(β1 )
(5) (6)
From its parameters R and β, we obtained relaxation times quantificationally. The KWW parameters for six kinds of ACF and the relaxation times are listed in Table 1. It shows that the end-to-end vector took the longest time, ∼960 ps in solution particularly. In the present study, our equilibrium time is longer than 1500 ps, and the sampling time is 3000 ps. Both of them are much larger the biggest τ measured here. Therefore, the duration length of the present MD run is long enough. The results from the previous two aspects show that the present simulation was moving in a normal way with reliable nonbonding interaction parameters and duration length. The following sections show our main contributions: to reveal the detailed solvation structures and to recognize the solvation mechanism for the polymer in solution. Detailed Solvation Structures of Biphenyl Molecule. In a solvation, specific interaction configurations must be formed between the solute and the solvent. Such structures are generally called solvation structures as was mentioned in the Introduction. In the present study, we approach the solvation structures in detail through two windows: the structure of biphenyl and the conformation of PE. We investigate the former in this section and the latter in the next section. The biphenyl molecules, which we focused on, are the biphenyls in the first solvation shell. The first solvation shell, as shown in other works,41-43 includes all the atoms on the PE chain and first neighbor atoms from the solvent, where the specific interaction configurations were expected to be found. The boundary of the first solvation shell was determined in terms of the relative distance of the biphenyl mass center and the
22722 J. Phys. Chem. B, Vol. 110, No. 45, 2006
Ji and Yang
Figure 4. Sketch map of classifying biphenyls in first solvation shell to three groups. Group I weighted 3%. Group II weighted 45%. Group III weighted 52%.
Figure 2. PCF for the mass center of biphenyl and carbon atom in PE chain.
Figure 5. Torsional distribution of two benzene rings in biphenyl. (0) Group I; (O) group II; and (4) group III. Dashed line for biphenyl in bulk. Figure 3. Scheme of two possible solvation structures.
nearest carbon atom in PE. It was investigated by pair correlation function (PCF). The PCF gives a measure of probability that there is an atom in a spherical shell with a distance r from the reference atom. In other words, this function measures a distance distribution between the solute and the solvent. It is generally expressed as
gAB(r) )
〈
V0NAB(r)
Vs(r)NABt
〉
(7)
v0
where NAB(r) is the number of AB pairs in the rth shell, NABt is the total number of AB pairs in the system, Vs is the volume of the rth shell, and V0 is the total volume of the system. The calculated PCF curve for the biphenyl mass center and the carbon atom in PE is shown in Figure 2. Clearly, there are three main regions in Figure 2 as labeled. Region 3 is about 9 Å, while region 1 is about 4.5 Å. Region 2 is about 6.7 Å. It is obvious that the second period region of region 1 should be located at ∼9 Å, where it is just region 3. This thus indicates that the first solvation shell includes regions 1 and 2. A valley appears around 8 Å, which implies that we can set the radius of the first solvation shell of the PE chain in biphenyl to 8 Å. The interaction configurations of the biphenyls in regions 1 and 2 probably come from the two models as shown in Figure 3 in the present study. That is to say, the biphenyl is approximately parallel with or vertical to the backbone of the PE chain. These suppositions will be confirmed next. To recognize the interaction configurations in detail, we classify all biphenyls in first solvation shell into three groups. The full statistical process of the classification is as follows. First, we measured the distance of the mass centers of each biphenyl to its nearest carbon atom in the PE chain. If the distance was smaller than 8 Å, the biphenyl was within the first solvation shell. Second, we divided all the biphenyls in the first solvation shell into three groups: group I has both end carbon atoms, which are through the long axis of the biphenyl molecule,
within 4.5 Å of the PE chain; group II has only one of the end carbon atoms within 4.5 Å of the PE chain; and group III has both end carbon atoms 4.5 Å from the PE chain. Their schemes are shown in Figure 4. We analyzed these three groups of interaction configurations in the first solvation shell. It was found that the population distribution of the three groups is that group I is 3%, group II is 45%, and group III is 52%. The detailed solvation structures of the biphenyl in the first solvation shell were considered in two ways: by the intramolecular structures and by the intermolecular structures. The dihedral angle between the two benzene rings in biphenyl is to be studied in the present work. Statistical distributions of the dihedral angles of the three groups and the bulk biphenyl are displayed in Figure 5. It is obvious that group II, group III, and the bulk biphenyl have similar distributions of the dihedral angle. Their peaks are overlapped, and the tops of the peaks are all present at ∼32°, while the distribution of group I is obviously different from the others. It is enhanced in the small angle side and reduced in the large angle side. Its peak is presented at ∼27°. This result shows a strong solvation effect of PE on biphenyls in group I: a close interaction leads to a 5° shift to the dihedral angle distribution of the biphenyls. This study shows that such tiny differences in the dihedral angle can be detected. But, on the rotational isomer state (RIS)15 approximation, the tiny difference can be omitted. The intermolecular interaction configuration of biphenyl, which we are interested in, is related to the orientation angle Φ between PE backbones and biphenyl principal axis. An illustration of the orientation angle is shown in Figure 6. To obtain the orientation angle, we measured the length of L and the distance difference
δm,n ) x(d2 - d1)2
(8)
Structures and Mechanism of Polymer Solvation
J. Phys. Chem. B, Vol. 110, No. 45, 2006 22723
Figure 6. Illustrations of angle Φ and δm,n.
Figure 8. Illustration of dyad selection.
TABLE 2: Conformation Fractions of PE Dyads in Various Casesa A I II III B
TT
TG
TG'
GT
GG
GG'
G'T
G'G
G'G'
0.418 0.452 0.398 0.313 0.398
0.131 0.121 0.142 0.152 0.124
0.123 0.111 0.129 0.142 0.125
0.112 0.130 0.134 0.145 0.127
0.040 0.022 0.036 0.047 0.040
0.008 0.002 0.008 0.012 0.010
0.123 0.139 0.120 0.141 0.126
0.009 0.001 0.006 0.008 0.010
0.035 0.020 0.028 0.041 0.039
a A and B for the whole chain in solution and in melt. I: Group I; II: group II; and III: group III.
Figure 7. Statistical distribution of the angle Φ for all biphenyls and each biphenyl group in first solvation shell.
Then, we obtained
φ ) arcsin(δm,n/L)
(9)
In terms of Φ dependence, the three groups have their own statistical distributions on the angle Φ as shown in Figure 7, as compared with all biphenyls in the first solvation shell. Besides abundance, there are significant differences between the four kinds of distributions. The triangle curve (group I) and square curve (group III) both have a top value near 0°, while there are two peaks existing in the circle curve (group II). They are about 11 and 61°, respectively. Through comparing the distribution of the three groups, one can find that the detailed interaction configurations in the first solvation shell possess the following orientation features: the biphenyls in group I prefer to be parallel to the chain backbone, those in group II prefer to have two special orientations with Φ around 11 and 61°, and those in group III have a broad distribution with a large population in a parallel orientation and a small one in a perpendicular orientation. On the whole, the overall feature of all biphenyls in the first solvation shell has been demonstrated by the solid line as shown in the Figure 7. It is obvious that except for the top value around 0°, there is a broad peak around 60°. This indicates that in the first solvation shell, there are two main orientation angles between the PE chain and the long axis of the biphenyl molecules. It is close to the suppositions related to Figure 3. Detailed Conformation Structures of Polymer Solute. In the present study, the detailed solvation structures of polymer solute were exposed to the local conformations of PE interacting with three group biphenyls. Related discussions here were organized by two points. One is to observe specific conformations caused by biphenyls in the three groups. The other is to test out differences in the conformational distribution of the end location from the middle location. As a general conformational classification,29 the torsion, φ, of a bond along the PE chain has three conformational states: T, G, and G′, and they are located at 120° < φ < 240°, 0°
