182
Langmuir 2003, 19, 182-185
Interactions and Assembly of Simple Solvated Polymers Models Laura J. Douglas Frink* and Andrew G. Salinger Sandia National Laboratories, Albuquerque, New Mexico 87185 Received July 10, 2002 In this paper we discusses the application of a three-dimensional implementation of a geometry-based nonlocal density functional theory to the interactions between rigid polymers in solution. Deoxyribonucleic acid (DNA) provides inspiration for the three simple model polymers (cylindrical, bead chain, and periodic) we consider; however the results are more generally applicable to rigid polymers such as liquid crystals. The calculations show that fine details of polymer surface structure play a critical role in determining the solvation energy landscape. This solvation energy landscape in turn controls routes for assembly of bundles of macromolecules. This approach provides a new level of molecular insight into interactions in solvated polymer systems with wide applications to both industrial and biological polymers.
1. Introduction Polymer solutions are critical to the processing of a wide range of materials (from thin films for optical applications1 to paper2). They are also critical to biological function as proteins are nearly always solvated in physiological conditions.3 In fact, recent simulatons of peptide folding have confirmed the importance of including explicit solvent in the calculations.4 Theories for solvated polymers usually begin from a solution point of view and focus on capturing polymer physics (polymer configurations, radius of gyration, etc.) while treating the solvent of the system with a single parameter that defines the strength of polymer-solvent interactions.5 In this paper we take an alternate approach that treats the polymers as surfaces that generate an external field in which the fluid molecules equilibrate. Polymers then interact through a solvent-mediated potential of mean force. We have calculated potentials of mean force between three model polymers with a novel three-dimensional (3D) nonlocal density functional theory code.6 This approach allows for detailed free energy calculations of the solvated interactions of locally rigid interacting polymers. The calculations presented here show how polymer geometry can affect the solvation energy landscape and, more specifically, how assembly of polymer bundles can occur via minor adjustments in the alignment of two polymer strands with respect to one another. 2. Theoretical Approach The statistical mechanics behind our approach has been detailed elsewhere.7,8 Briefly, the semigrand ensemble for a solution of N macromolecules in a solvent with known chemical potential, µ, at a temperature, T, and volume, V, is (1) Lee, J. W.; Wang, C. S.; Price, G. E.; Husband, D. M. Polymer 1997, 38, 1403. (2) Nada, A. M. A.; Abdelhakim, A. A.; Mohamed, E. S.; Badran, A. S. J. Elastomers Plast. 1998, 30, 29. (3) Cheng, Y. K.; Rossky, P. J. Nature 1998, 392, 696. (4) Daura, X.; Mark, A. E.; van Gunsteren, W. F. Comput. Phys. Commun. 1999, 123, 97. (5) de Gennes, P.-G. Scaling Concepts in Polymer Physics; Cornell University Press: Ithaca, NY, 1979. (6) Frink, L. J. D.; Salinger, A. G. J. Comput. Phys. 2000, 159, 407. (7) Douglas, L. J.; Lupkowski, M.; Dodd, T. L.; van Swol, F. Langmuir 1993, 9, 1442. (8) Frink, L. J. D.; van Swol, F. J. Chem. Phys. 1994, 100, 9106.
ZNVµT )
{
∞
zn
∫V dRN e-βU (R ) ∑ n! ∫V drn e-βU (R ,r ) n)0 N
N
n
N
n
}
(1)
where z ) eβµ/Λ3 is the solvent activity, β ) 1/kBT, kB is the Boltzmann constant, and Λ is the de Broglie wavelength. Integrals are taken over all possible configurations of the N surfaces (denoted RN) and the n solvent molecules (denoted rn). The potential energy of direct surface-surface interactions is given by UN while the energies of both solventsolvent and solvent-surface interactions are summed in Un. Identifying the inner integral in eq 1 as the grand partition function, ΞµVT, of a fluid in an external field of N fixed surfaces, the mixture partition function may then be written in terms of the grand potential
Ω(RN;µ,T) ) -β-1 ln ΞµVT
(2)
Substituting for ΞµVT, eq 1 becomes
ZNVµT )
∫V dRN e-βU (R ) e-βΩ(R ;µ,T) N
N
N
(3)
Taking two hard rigid macromolecules (UN ) ∞ if overlapping, otherwise, UN ) 0), the potential of mean force (PMF) acting between the macromolecules is
W(R) ) Ω(R;µ,T) - Ω(∞;µ,T)
(4)
in the limit of two surfaces separated by R in an otherwise infinitely dilute solution of the macromolecules (i.e., FM f 0). Providing that the Gibbs dividing surface is chosen to be identical at all surface separations, the PMF may also be written in terms of the surface free energy, W(R) ) Ωs(R;µ,T) - Ωs(∞;µ,T), where Ωs ) Ω(R;{F(r)},T) Ω(R;{Fb},T), where Fb is the bulk fluid density associated with the known µ and F(r) is the equilibrium nonuniform density distribution of fluid particles in the external field of the N macromolecules. The surface free energy is directly related to the solventmediated (or solvation) force acting on the polymer strands.
10.1021/la026212m CCC: $25.00 © 2003 American Chemical Society Published on Web 12/07/2002
Solvated Polymer Models
Langmuir, Vol. 19, No. 1, 2003 183
Given hard impenetrable surfaces, this solvation force may be calculated from the sum rule9
fx ) -
∂Ωs ) ∂R
∫ F(rs)nx drs
(5)
where ∫ drs indicates a integral over the surface of the macromolecule and n is the unit normal to the surface. The inhomogeneous solvent density distribution, F(r), and the surface free energies, Ωs, and consequently the potential of mean force, W, can be calculated with molecular simulation or density functional theory (DFT). DFTs are based on the functional minimization of the grand free energy, Ω[F(r)] with respect to the density distributions, F(r) at constant temperature, T, and fluid chemical potential, µ.
( ) δΩ δF(r)
)0
(6)
T,µ
The particular functional formulation we apply here is very accurate for inhomogeneous hard sphere fluids and was originally developed by Rosenfeld.10 While most DFT calculations have considered geometries with considerable symmetry (slit-pores, cylindrical pores, spherical cavities, etc.); the polymer calculations here will require a full 3D DFT solution. Go¨tzelmann et al. have shown that when an accurate equation of state is known for the macromolecule-solvent mixture, PMFs can be found by taking the explicit limit of FM f 0.11 While their approach shows great promise for certain geometries (e.g., spherical colloids near flat surfaces), it cannot be easily extended to the polymer-solvent mixtures considered here because reliable free energy functionals for the binary polymersolvent mixtures are not available. In addition while the sphere-plate geometry has the advantage of one surface with simple symmetries (the plate), which results in a 1D calculation of solvent in its fixed external field, there are no simple surfaces in the interaction of two polymer strands. Treating one polymer strand as the fixed surface would still require a 3D calculation. Our three-dimensional numerical implementation is based on a Newton’s method solve of the system of equations, and convergence is usually obtained in less than 20 Newton iterations.6 Each of the 3D polymer solutions presented here was obtained in approximately 3 min of CPU time on 50 processors of the ASCI-Red (Intel pentium 333 MHz chips) computer at Sandia National Laboratories.
Figure 1. Two interacting DNA strands with atomic coordinates taken from the protein data bank.
We restrict the current discussion to hard sphere fluids in contact with hard rigid polymers. While other physical effects such as van der Waals, polarization, or Coulomb forces may be critical in many systems, we consider only volume exclusion in order to demonstrate our approach, the benefit of having a 3D numerical method for arbitrary surfaces, and the richness of the physics that arises even from volume exclusion alone. The stiffness of the polymer models discussed here is not a realistic representation of all polymers; however there are some important exceptions. One example is deoxyribonucleic acid (DNA), which is rigid over long persistence lengths, λp, due to the double helix that
composes its backbone. Due to this rigidity, DNA is in fact often treated as a rigid cylindrical polymer in molecular investigations.12-14 More generally, every polymer is rigid on length scales smaller than its characteristic persistence length.5 Thus part of our discussion here involves identifying the λp associated with strong (>kT) solvation energy barriers. Using DNA as a trial system (see Figure 1), we consider three possible coarse models based on the overall topology of the chains. They include cylindrical, bead chain, and periodic models. Sketches of the three models can be found in Figure 2. Taking the z axis down the long axis of a polymer chain, the cylindrical polymer is described by Rp(z) ) 1.5σ, where Rp is the radius of the polymer chain and σ is the diameter of the solvent particle. The bead chain polymer is composed of spheres with radius Rp ) 1.5σ, where the centers of the spheres along the z axis are separated by 3σ. The periodic polymer has a radius Rp(z) ) ro + A cos(2πz/λ), where the reference radius ro ) 0.75σ, the amplitude A ) 0.75σ, and the period λ ) 3σ. While all three polymer models exhibit surface curvature perpendicular to the long axis of the polymer, the bead chain and periodic polymers also exhibit structure parallel to the polymer axes. Thus, we have calculated potentials of mean force as a function of both relative alignment and surface separation for these two cases. Potentials of mean force for the three solvated polymer models are shown in Figure 2. In all cases, solvent packing leads to oscillatory potentials as a function of the interaxial separation of the polymers. The magnitude of the oscillations is largest for the parallel cylinders in Figure 2A where densities are uniform along the long axis of the polymer strands. However, in all cases, the free energy peaks are substantial in comparison with kT. Consider the peak at R/σ ) 3.6 in the case of the aligned bead-chain polymer. It’s magnitude per unit length is Wσ/L ≈ 0.25kT. The persistence length needed for this energy barrier to be W ≈ kT is λp ≈ 4σ (or λp ≈ 2 nm for a σ ) 0.5 nm solvent
(9) D. Henderson, Ed. Fundamentals of Inhomogeneous Fluids; Marcel Dekker: New York, 1992. (10) Rosenfeld, Y. Phys. Rev. Lett. 1989, 63, 980. (11) Go¨tzelmann, B.; Roth, R.; Dietrich, S.; Dijkstra, M.; Evans, R. Europhys. Lett. 1999, 47, 398.
(12) Lyubartsev, A. P.; Nordenskiold, L. J. Phys. Chem. 1995, 99, 10373. (13) Gro¨bech; Jensen, N. Mashl, R. J.; Bruinsma, R. F.; Gelbart, W. M. Phys. Rev. Lett. 1997, 78, 2477. (14) Allison, S. A.; Mazur, S. Biopolymers 1998, 46 359.
3. Discussion of Results
184
Langmuir, Vol. 19, No. 1, 2003
Frink and Salinger
Figure 3. The density distribution as a function of position (x, y in units of σ) in a slice of constant z from a three-dimensional calculation on the misaligned bead-chain model (see sketch in Figure 2B). The z slice was taken at the point where two beads on the left chain come together. The white regions include both the polymer volume and solvent exclusion zones due to the hard interactions between the polymer and solvent particles. The maximum density in the three-dimensional solution, Fσ3 ) 10.7, is found in the plane of this figure at x/σ ) 6.375, y/σ ) 4.5, 6.5.
Figure 2. The potential of mean force as a function of the interaxial separation, R, of two cylindrical (A), bead-chain (B), and periodic polymers (C). In B and C, the solid lines show the cases where the polymers are perfectly aligned, while the dashed lines show cases where the polymers are misaligned as shown in the sketches. The bulk fluid density away from the polymers was Fbσ3 ) 0.63 for all cases.
molecule). Since this requisite persistence length is quite short, these results should be applicable to many polymer systems. The bead-chain case in Figure 2B demonstrates that the details of the surface structure of the polymer play an important role in assembly of the strands. Consider the energy barriers experienced by two polymer strands as they come together from infinite separation. If the strands remain in a perfectly aligned state (solid curve in Figure 2B), they will need to overcome the substantial solvation barriers mentioned above. However by shifting their relative alignment (dashed curve in Figure 2B) at the appropriate separations, the solvation barriers can be avoided altogether since the maxima in the aligned state very nearly correspond to the minima for the unaligned state. As a corallary, the forces experienced by the beadchain polymers may be always attractive from infinite separation to contact provided that the polymers can freely change their relative alignment. The only way for the cylindrical polymers to avoid the large solvation barriers would be for the cylinders to rotate away from the parallel position. Such a rotation would be kinetically unfavorable in comparison with the small adjustments needed when surface structure is present. While the surface structure on the bead-chain polymers provides a low-energy route to assembly by avoiding solvation energy barriers, this is not the case for the periodic polymer as shown in Figure 2C. Here, the unaligned and aligned cases have minima and maxima at similar separations. Nevertheless, the solvation energy
landscape is significantly affected by the surface structure on the periodic polymers. Specifically, the period of solvation oscillations increases from 1σ (Figure 2A,B) to approximately 1.5σ (Figure 2C). In addition, the global free energy minimum is now found at R/σ ) 2.75, where there is a layer of fluid between the polymer strands. This solvated assembly predicted for the periodic polymers is in contrast to both cylindrical and bead-chain cases (Figure 2A,B) where the global free energy minimum is found at contact. In these two cases, osmotic exclusion of the hard-sphere solvent results in the strong attractive forces at contact.11 Thus while solutions composed of any of the three polymers would self-assemble in to tightly ordered arrays, the periodic polymer assembly will be the most highly solvated. Predicting the existence of solvated assemblies from the theory is important for further studies of swelling of ordered phases and investigations of the forces that stabilize highly solvated crystals. For example, swollen DNA assemblies have been observed in experiments,15,16 and protein crystals may contain 50% water by volume.17,18 For all of the cases in Figure 2 except the unaligned bead-chain polymer, the potential curves are drawn to the point of closest approach of the polymer strands. However, it was not possible to obtain solutions for the case of the unaligned bead-chain polymer when R/σ e 3.5. The physical origin of the numerical instability is local crystallization (or physical binding) of solvent in the annuli around the points where two beads on one chain meet (see Figure 3). In this case, narrow steep density peaks must be resolved. This becomes impossible for a direct solve on a fixed mesh grid using linear basis functions to represent the density profiles. This problem may be exaserbated in our calculations by using Rosenfeld’s original functional that is known to exhibit incorrect zero-dimensional crossover behavior.19 Due to the numerical limitations of the collocation approach, a Gaussian basis approach used (15) Rau, D. C.; Parsegian, V. A. Biophys. J. 1992, 61, 260. (16) Frink, L. J. D.; van Swol, F. Colloids Surf., A 2000, 162, 25. (17) Smith, P. E.; Pettitt, B. M. J. Phys. Chem. 1994, 98, 9700. (18) Schoenborn, B. P. J. Mol. Biol. 1988, 201, 741. (19) Rosenfeld, Y.; Schmidt, M.; Lo¨wen, H.; Tarazona, P. J. Phys.: Condens. Matter 1996, 8, L577.
Solvated Polymer Models
in DFT calculations of crystals20 may be needed for cases where strong solvent localization is important. 4. Summary To summarize, this paper presents detailed threedimensional density functional theory calculations of the solvation energy landscape between interacting solvated polymers. This approach is quite different from most treatments of polymer solutions where the solvent is a continuum. For each of the three polymer surface topologies considered, the theory predicts that solvent packing can be an important effect even for polymers with rather small (e.g., λp ≈ 2 nm) persistence lengths. Furthermore, it was shown how the fine details of surface topology may affect the solvation energy landscape through which the polymers interact. This solvation energy landscape in turn (20) Ohnesorge, R.; Lo¨wen, H.; Wagner, H. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 1994, 50, 4801.
Langmuir, Vol. 19, No. 1, 2003 185
controls both the available pathways to assembly of bundled arrays of polymer strands and the degree of solvation in the self-assembled state that ultimately forms. While the link between polymer topology, solvent properties, and assembly is clearly complex, these calculations show that it is possible to assess these fine molecular details using a DFT-based approach. Our current work involves replacing the rather crude polymer models discussed here with models based on known atomic configurations of complicated polymers such as DNA and proteins. Acknowledgment. The authors thank Mark Stevens and Jeff Weinhold for many helpful discussions. Sandia is a multiprogram laboratory operated by Sandia Corporation, a Lockheed Martin Company, for the United States Department of Energy under Contract DEAC04-94AL85000. LA026212M
