- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Analysis of interactions between finite-sized particles...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Analysis of interactions between finite-sized particles and terminally attached polymer using numerical.. 1997
pdf
Page Metadata
Item Metadata
Title | Analysis of interactions between finite-sized particles and terminally attached polymer using numerical self-consistent-field theory |
Creator |
Steels, Bradley Michael |
Date Created | 2009-04-28 |
Date Issued | 2009-04-28 |
Date | 1997 |
Description | Non-specific protein adsorption is known to influence many processes. For example, fouling of food processing equipment or fouling of biomedical implants often occurs because of strong protein adsorption. Protein adsorption on these and other surfaces, such as biomaterials, can however be reduced by attaching polymer chains by one end to the surface. Terminally attached polymer chains have also found application in size-based chromatographic separations. The goal of this thesis is to improve our understanding of the way in which solute molecules interact with terminally attached polymer chains, and to use this knowledge to predict optimum system conditions for minimizing protein adsorption. We develop a numerical selfconsistent- field lattice model, based on an earlier model of Scheutjens and Fleer [1979; 1980], to calculate theoretical results for the polymer density distribution of isolated and interacting chains around a solute particle positioned at a fixed distance from a surface. In addition, the excess energy required to move the particle into the polymer chains (interaction energy) is calculated using a statistical mechanical treatment of the lattice model. The effect of system variables such as particle size, chain length, surface density and interaction parameters on density distributions and interaction energies is also studied. The model is first applied to the compression of a single polymer chain by a disc of finite radius. Results are compared to predictions of a previous scaling thermodynamic model by Subramanian et al. [1995]. We see no first order partial-escape transition as reported by Subramanian et al. Instead, we find that the energy required to compress the chain increases monotonically, becoming independent of chain length at very close compression Calculations for the interaction of a solute particle with a surface covered by many polymer chains (a brush) shows that the polymer segments will fill in behind the particle quite rapidly as it moves toward the surface. When there is no strong energetic attraction between the polymer and solute we predict that the interaction energy will be purely repulsive upon compression due to losses in conformational entropy of the polymer chains. Above a critical chain length, which depends upon particle size, a maximum in the force required to move the particle toward the surface is observed due to an engulfment of the particle as chains attempt to access the free volume behind the particle. If an attraction exists between the polymer and solute, such that a minimum in the interaction energy is seen, the optimum conditions for solute repulsion occur at the highest surface density attainable. Long chain length can lead to increased solute concentration within the polymer layer due to the fact that more favourable polymer-solute contacts are able to occur than with short chains at a similar entropic penalty. |
Extent | 6736785 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
File Format | application/pdf |
Language | Eng |
Collection |
Retrospective Theses and Dissertations, 1919-2007 |
Series | UBC Retrospective Theses Digitization Project [http://www.library.ubc.ca/archives/retro_theses/] |
Date Available | 2009-04-28 |
DOI | 10.14288/1.0058581 |
Degree |
Master of Applied Science - MASc |
Program |
Chemical and Biological Engineering |
Affiliation |
Applied Science, Faculty of |
Degree Grantor | University of British Columbia |
Graduation Date | 1998-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
URI | http://hdl.handle.net/2429/7658 |
Aggregated Source Repository | DSpace |
Digital Resource Original Record | https://open.library.ubc.ca/collections/831/items/1.0058581/source |
Download
- Media
- ubc_1998-0056.pdf [ 6.42MB ]
- Metadata
- JSON: 1.0058581.json
- JSON-LD: 1.0058581+ld.json
- RDF/XML (Pretty): 1.0058581.xml
- RDF/JSON: 1.0058581+rdf.json
- Turtle: 1.0058581+rdf-turtle.txt
- N-Triples: 1.0058581+rdf-ntriples.txt
- Citation
- 1.0058581.ris
Full Text
Analysis of Interactions between Finite-Sized Particles and Terminally Attached Polymer using Numerical Self-Consistent-Field Theory by BRADLEY MICHAEL STEELS B.Eng., Royal Military College of Canada, 1995 , A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTERS OF APPLIED SCIENCE in THE FACULTY OF GRADUATE STUDIES Department of Chemical and Bio-Resource Engineering and The Biotechnology Laboratory We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA December 1997 © Bradley Michael Steels, 1997 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of CJA&M iCW OA^A BlOfCSourC^ LZ<^-C\ * < * - * £ > * A - ' ^ The University of British Columbia Vancouver, Canada DE-6 (2/88) ABSTRACT Non-specific protein adsorption is known to influence many processes. For example, fouling of food processing equipment or fouling of biomedical implants often occurs because of strong protein adsorption. Protein adsorption on these and other surfaces, such as biomaterials, can however be reduced by attaching polymer chains by one end to the surface. Terminally attached polymer chains have also found application in size-based chromatographic separations. The goal of this thesis is to improve our understanding of the way in which solute molecules interact with terminally attached polymer chains, and to use this knowledge to predict optimum system conditions for minimizing protein adsorption. We develop a numerical self- consistent-field lattice model, based on an earlier model of Scheutjens and Fleer [1979; 1980], to calculate theoretical results for the polymer density distribution of isolated and interacting chains around a solute particle positioned at a fixed distance from a surface. In addition, the excess energy required to move the particle into the polymer chains (interaction energy) is calculated using a statistical mechanical treatment of the lattice model. The effect of system variables such as particle size, chain length, surface density and interaction parameters on density distributions and interaction energies is also studied. The model is first applied to the compression of a single polymer chain by a disc of finite radius. Results are compared to predictions of a previous scaling thermodynamic model by Subramanian et al. [1995]. We see no first order partial-escape transition as reported by Subramanian et al. Instead, we find that the energy required to compress the chain increases monotonically, becoming independent of chain length at very close compression. ii Calculations for the interaction of a solute particle with a surface covered by many polymer chains (a brush) shows that the polymer segments will fill in behind the particle quite rapidly as it moves toward the surface. When there is no strong energetic attraction between the polymer and solute we predict that the interaction energy will be purely repulsive upon compression due to losses in conformational entropy of the polymer chains. Above a critical chain length, which depends upon particle size, a maximum in the force required to move the particle toward the surface is observed due to an engulfment of the particle as chains attempt to access the free volume behind the particle. If an attraction exists between the polymer and solute, such that a minimum in the interaction energy is seen, the optimum conditions for solute repulsion occur at the highest surface density attainable. Long chain length can lead to increased solute concentration within the polymer layer due to the fact that more favourable polymer-solute contacts are able to occur than with short chains at a similar entropic penalty. iii Table of Contents ABSTRACT II TABLE OF CONTENTS IV TABLE OF FIGURES VI LIST OF TABLES XI ACKNOWLEDGEMENTS XII 1 INTRODUCTION 1 1.1 Overview 1 1.2 Non-Specific Protein-Surface Interactions 4 1.3 Thesis Objective 6 1.4 Terminally Attached Polymers at Interfaces 7 1.4.1 Distributions of Terminally Attached Chains 10 1.4.1.1 Experimental Measurements and Scaling Predictions 10 1.4.1.2 Advanced Brush Models 16 1.4.1.3 Brushes with Free Polymer in Solution 19 1.4.1.4 Brush Distributions - Summary 20 1.4.2 Compression of Terminally Attached Polymers 22 1.4.3 Particle Interactions with Terminally Attached Polymers 28 1.4.3.1 Experimental Studies 28 1.4.3.2 Modeling of Brush-Particle Interactions 30 2 MODEL OUTLINE 33 2.1 Origins of Nonideal Mixing in Liquid Solutions 33 2.1.1 Entropy of Mixing Polymer in Solvent 34 2.1.2 Enthalpy and Gibbs Energy of Mixing 36 2.2 Self-Consistent Mean-Field Theory of Scheutjens and Fleer 38 2.2.1 Cylindrical Lattice Geometry 41 2.2.2 Combinatorial Entropy 43 2.2.3 Internal Energy 45 2.2.4 Equilibrium Segment Density Distribution 46 iv 2.2.5 Chain Grafting 48 2.2.6 Interaction Energy 51 2.2.7 Numerical Solution 52 2.3 Evaluation of Model Parameters 57 2.3.1 Scaling the Lattice Model 57 2.3.1.1 Lattice Unit and Chain Length 58 2.3.1.2 Polymer Surface Concentration 61 2.3.1.3 Particle Size 62 2.3.2 Measuring Interaction Parameters 64 2.3.2.1 Differential-Vapour-Pressure and Osmometry 65 2.3.2.2 Low-Angle Laser-Light-Scattering (LALLS) 66 3 ISOLATED TERMINALLY-ATTACHED CHAINS 71 3.1 Motivation 71 3.2 Scaling Theory of Mushroom Compression 72 3.3 Analysis of Mushroom Compression using Numerical SCF Theory 73 3.4 Summary 83 4 BRUSHES AND BRUSH-PARTICLE INTERACTIONS 84 4.1 Outline of Calculations 84 4.2 Brush Density Distributions 85 4.3 Brush-Particle Interactions 89 4.3.1 2D Brush Density Distributions 89 4.3.2 Brush-Particle Interaction Energy 94 4.3.3 Surface Mobile Brushes 101 4.3.4 Brush-Particle Attraction 104 NOMENCLATURE 113 REFERENCES 117 APPENDIX A 127 v Table of Figures Figure 1.1 Several possible polymer density profiles are shown for terminally attached polymer chains, (a) The Alexander-de Gennes brush assumes a step profile, (b) Analytical SCF theory predicts a parabolic brush, (c) Real terminally attached polymer will exhibit a maximum slightly away from the surface, which advanced numerical models also predict, (d) A strongly adsorbing surface can lead to a pancake structure 14 Figure 2.1: Comparison between ideal entropy of mixing and combinatorial entropy of mixing given by Flory-Huggins theory for a polymer chain of 1000 segments 36 Figure 2.2: Cartesian (a) and cylindrical (b) lattice geometries used in the self-consistent-field model. Areas treated with mean-field averaging are shown in grey as well as the grafting surface in (b) 39 Figure 2.3: The probability that segment s ends in radial (r,z) is found by combining the probability of two shorter chains ending in that radial. This involves combining two step- weighted walks, one for 1,2,...s segments and one for N,N-1,...s 49 Figure 2.4: Non-linear curve fit of data taken from Tyn and Gusek [1990] for use scaling particle size to protein molecular weight 64 Figure 3.1 Schematic diagram of a single end-grafted chain in (a) uncompressed (b) compressed, and (c) escaped states. The chain grafting point indicated by • is at lattice position (z=0,r=0). Characteristic dimensions include: (i) disc position, Z^ c , (ii) disc radius, Rdisc (iii) r.m.s. radius of segments from cylinder axis, Rrms and, (iv) Flory radius, RF. 71 Figure 3.2 Fractional segment distribution of a single end-grafted chain in the radial (r) (dash- dot line) and normal (z) (solid line) directions with iV=100 segments. Dashed line shows the root mean squared (rms) distance of segments from the center axis, Rrms, and dotted line shows the Flory radius, RF 75 vi Figure 3.3 Segment density distribution of a single end-grafted chain with N=100. The chain is shown in (a) with no disc and in (b) with a disc of radius, Rdisc=20, sitting in layer 2 compressing the mushroom 76 Figure 3.4 Root mean squared distance of segments from the center of the cylinder, Rtms, is shown as the disc compresses a mushroom with N=100. Rdisc=20 (solid), Rdisc=30 (dashed) 77 Figure 3.5 Fractional segment distribution in the radial direction as a mushroom with N=75 is compressed by a disc with Rdisc=20. The curves represent the disc in positions of z=16, 8, 5, 4, 3, and 2 78 Figure 3.6 Root mean squared distance of segments from the center of the cylinder, Rims, is shown as the disc of RdiSc=20 compresses a mushroom with N=75 79 Figure 3.7 Energy required to compress a mushroom with a disc of radius, Rdisc=20, for several chain lengths: N=20 (dashed), N=50 (dotted), N=100 (dot-dash), N=200 (solid) 80 Figure 3.8 Fraction of segments covered (dotted) and escaped (solid) from under a stationary disc of radius 20, extending from z=3 to the end of the cylinder (Zmax), as chain length is varied. 81 Figure 3.9 Chemical potential of a single end-grafted chain under a stationary disc of radius 20, that extends from z=3 to the end of the cylinder ( Z m a x ) , as a function of the number of grafted chains 82 Figure 3.10 The change in chemical potential with respect to the number of grafted chains (at nchain = 1) as a function of chain length. This value always remains positive indicating stability against phase separation 82 Figure 4.1 Schematic diagram showing a solute molecule interacting with polymer chains end- grafted to a surface. The centerline represents the axis of a cylindrical coordinate system, and the rings on the surface represent spatial increments in the radial direction 84 vii Figure 4.2 Brush density distributions for various polymer-solvent interaction parameters. N = 50 segments, o = 5%, = jtf>s = 0 85 Figure 4.3 Brush density profiles for various polymer-surface interactions are shown. N = 50 segments, cr = 5%, %os = Xbo = 0 87 Figure 4.4 Segment density profiles for brushes of different lengths. Surface density is 10% and the polymer-solvent interaction parameter %b0 1S 0.5. All other interactions are athermal 88 Figure 4.5 Brush density profiles for various surface densities. Chain length is 100 segments and^bo = 0.5 88 Figure 4.6 Polymer segment density distributions in the cylindrical lattice are indexed in both the z and r directions. The chains of length N=50 with a grafting density of c=0.1 are interacting with a particle of dimensions Rp=3 and Lp=3 (3 by 3). Interaction parameters are set as Xbo=0.4 and Xpo=0.4 (unspecified interaction parameters are always zero). The density distribution is shown for a particle position (a) outside the brush, (b) at layer Zp = 10, (c) zp = 7, and (d) Zp = 5 90 Figure 4.7 The 2-dimensional segment density distributions for four different chain lengths are shown around a 3 by 3 particle placed with the leading edge at Zp = 5. Interaction parameters are specified as %bo=0.4 and Xpo-0.4. Chain lengths are (a) N - 15, (b) N = 25, (c) N = 40 and (d) N = 60 segments 92 Figure 4.8 The brush distribution of N=50 segment chains is shown around a 5 by 5 particle at four different grafting densities. Interaction parameters are Xbo=0.4 and Xpo=0.4. Grafting densities are (a) a=0.05, (b) a=0.1, (c) a=0.2 and (d) a=0.3 93 Figure 4.9 Brush-particle interaction energy curves are shown for chains of N=50 with a graft density of a=0.1 interacting with various sized particles. Interaction parameters are specified as Xbo=0.4 and xPo=0.4. Particle size has been varied from Rp=l, Lp=l up to Rp=7, Lp=7. The dashed line shows the inflection point of each curve that shows one 94 viii Figure 4.10 Interaction energy curves are shown for a particle of Rp=3, Lp=3 interacting with different length polymer chains of constant grafting density a - 0.1. Interaction parameters are Xbo-0.4 and xPo=0.4. The dashed line shows the inflection point of each curve 95 Figure 4.11 The critical chain length Ncmcai for which an inflection point is observed is shown as a function of Rp for the system with a = 0.1, the previous interaction parameters, and particles with a cylindrical geometry where LP = RP 96 Figure 4.12 F and the rate of segment displacement from the total volume the particle moves through are shown as a function of zp for a 3 by 3 particle moving into a brush with N = 50 and cr=0.1 98 Figure 4.13 Fmax is shown as a function of particle size for a brush with N = 50 and a 10% graft density. The maximum force appears to scale closely with particle volume 100 Figure 4.14 Interaction energy curves are shown for a particle of Rp=3, Lp=3 interacting with 50 segment polymer chains of varied grafting density. Interaction parameters are Xbo=0.4 and XPo=0.4. The dashed line shows the inflection point of each curve 101 Figure 4.15 Brush distribution around a 3 by 3 particle is shown for surface non-mobile and mobile brushes. Brushes consist of N=40 segment chains at 10% surface density and interaction parameters used were Xbo~0.4 and xPo=0.4. The distributions for non-mobile chains around a particle in (a) layer zp = 4 and (b) zp = 2 are shown; distributions of mobile chains around a particle in (c) zp = 4 and (d) zp = 2 are also shown 102 Figure 4.16 Interaction energy curves are shown for a 3 by 3 particle interacting with 40 segment polymer chains that are surface mobile (dashed) and non-mobile (solid). Interaction parameters are Xbo=0.4 and xPo=0.4 103 Figure 4.17 Brush-particle interaction energy curves are shown for five different sized particles interacting with chains of length N=50 and grafting density a=0.1. Interaction parameters are Xbo=0.5, xPo=0.5 and Xb P = -0.5 105 ix Figure 4.18 Brush-particle interaction energy curves are shown for a 5 by 5 particle penetrating chains of length N=50 at four different grafting densities. Interaction parameters are Xbo=0.5, xPo=0.5 and Xb P= -0.5 106 Figure 4.19 Brush-particle interaction energy curves are shown for a 3 by 3 particle penetrating 1 brushes of varied chain length, each at 10% grafting density. Interaction parameters are Xbo=0.5, xPo=0.5 and Xb P= -0.5 108 Figure 4.20 Layer partition coefficient Kp(zp) for various particle sizes interacting with 50 segment chains at 10% surface density. Interaction parameters are Xbo=0.5, xPo=0.5 and Xb P= -0.5. 110 List of Tables Table 1.1 Several polymers with their monomers and repeat units 8 Table 1.2 Schematic of different copolymer structures 8 Table 2.1: Experimental properties that require direct comparison with lattice values 58 Table 2.2: Polymer-solvent interaction parameters for PEG and Dextran in water at 25 °C calculated from light-scattering data of Haynes et al. [1993] 69 Table 2.3: Protein-solvent interaction parameters calculated from light-scattering data for bovine serum albumin (BSA) in water at 25 °C 69 Table 2.4: Protein-solvent interaction parameters calculated from light-scattering data for lysozyme in water at 25 °C 69 Table 2.5: Protein-solvent interaction parameters calculated from light-scattering data for ct- chymotrypsin in water at 25 °C 70 Table 2.6: Protein-polymer interaction parameters X23 calculated from light-scattering data in potassium phosphate buffer at 25 °C 70 xi Acknowledgements Above all, I would like to thank my thesis supervisor Dr. Charles Haynes, for offering me this very interesting project. His enthusiasm and expertise have made this work possible. I am very grateful to Dr. Don Brooks and Dr. Bruce Bowen, whose courses in biophysical chemistry and numerical programming methods provided me a framework upon which this work was built. I would also like to thank Dr. Frans Leermakers for his expert help with the modeling and I am indebted to Jurgen Koska for his help with computer code and math. I must also thank all of the people that I have had the pleasure of knowing and working with in the Biotechnology Laboratory at UBC, especially those in the Bioprocess Engineering Group. This very outgoing group of people has been inspirational and supportive. I am thankful for the many wonderful friendships I have made during this work. I am also very grateful to NSERC and the Canadian Forces, whose support made this work possible. xii 1 Introduction 1.1 Overview Properties of a solution at an interface are different from those in bulk solution. This is especially true for solutions containing globular proteins, since the solubility and solution structure of proteins is determined by a delicate balance of intermolecular and intramolecular interactions. Introduction of an interface can disrupt this balance of forces, leading to adsorption and the possibility of a change in protein conformation. Protein adsorption is a complex process, which remains only partially understood [Haynes and Norde, 1994; Andrade, 1986]. Increasing our knowledge of protein-surface interactions is motivated by a desire to understand how these interactions may be controlled. Therapeutically effective contact of foreign surfaces with the human biosystem usually requires careful design of the sorbent properties. In some applications, it may be desirable for proteins to interact strongly with a surface, while in others adsorption must be moderate and controllable. Often, when a material is engineered for use in a biological system, protein adsorption is undesirable. When a foreign material contacts blood, plasma proteins such as serum albumin and fibrinogen form an adsorbed layer through a process called opsonization. It has been shown that nonspecific adsorption of proteins can lead to structural changes in the protein molecules [Haynes and Norde, 1994]. It is believed that these changes in conformation of serum proteins then trigger proteolytic cascades or activate platelet receptors, which leads to thrombus formation [Elbert and Hubbel, 1996; Bamford and Al-Lamee, 1992]. This has practical implications in design of all biomaterials, such as dialysis tubing, vascular grafts or artificial heart valves, since thrombus formation can lead to occlusion or embolus formation. 1 Protein adsorption in biosystems also has implications in site-directed drug delivery or gene therapy [Davis, 1997]. Adsorption of plasma proteins on the surface of a foreign material may lead to targeting of the material by the immune system, or more specifically by the mononuclear phagocytic system (MPS) [Woodle and Lasic, 1992]. Drug delivery vessels, such as liposomes or polymer microspheres, may be rapidly cleared from circulation without having time to reach the desired target within the host. This is especially significant in cases where a therapeutic is designed for second or third-order targeting, where it must reach a specific cell type or specific location within a cell, respectively. Unwanted protein adsorption is also observed in more general situations. Plaque formation on teeth and fouling of surfaces such as contact lenses, food processing equipment and naval equipment results from strong protein-surface interactions. The nonspecific adsorption of protein can also interfere with results in high-sensitivity bioassays in which fluids contact glass or polymeric surfaces [Malmsten and Van Alstine, 1996]. In some cases, protein adsorption can be a desired event. Because of their amphiphilic nature, proteins may be used as surfactant-like molecules to stabilize creams or emulsions at a liquid-liquid interface. Separation or purification of proteins may also be carried out using controlled protein-surface interactions. Hydrophobic interaction chromatography (HIC) makes use of hydrophobic interactions between protein surfaces and a solid support matrix. Under the correct solution conditions it is possible to adsorb target proteins preferentially from a mixture and then elute them under different solution conditions (i.e. lower salt concentration). Similarly, ion exchange chromatography makes use of electrostatic interactions between proteins and the support matrix. Protein purification processes are very important since therapeutic proteins are 2 often produced in biological fermentations, where the product concentration is typically quite low and the number and levels of impurities relatively high. A common step in protein purification is precipitation, which may be initiated by addition of water-soluble polymer to a protein-containing solution. Polyethylene glycol (PEG), dextran and polyvinyl alcohol (PVA) are often used as precipitating agents because they do not interact strongly with most proteins and are approved for biological applications by the FDA. The weak interaction of these water-soluble polymers with globular proteins has motivated their application in surface modifications for improved biocompatibility. Attaching PEG or dextran to synthetic surfaces has been shown to reduce significantly the adsorption of proteins and the adhesion of platelets and cells [Schroen et al, 1995; Llanos and Sefton, 1993: Merril, 1992; Bergstrom et al., 1992; Gombotz et al, 1992; Desai and Hubbel, 1991]. Presentation of PEG on the surfaces of polystyrene microspheres or liposomes has also been shown to reduce clearance of the particles from circulation [Dunn et al, 1994; Woodle and Lasic, 1992]. In both cases, the susceptibility of a surface to protein fouling can be controlled by the introduction of a water- soluble polymer to the interface. A new generation of size-exclusion chromatography, which we call entropic interaction chromatography, is based on end-grafting of water-soluble polymer chains to a solid matrix. Traditional size exclusion chromatography makes use of macroporous particles into whose pores some solute molecules may diffuse while others are excluded. Thus, the effective column volume is smaller for large molecules and they elute more quickly. A similar size-based separation is achieved in entropic interaction chromatography [Brooks and Muller, 1996]. In this case, smaller solute molecules sample a larger column volume as a result of lower steric repulsion by the polymer chains. 3 1.2 Non-Specific Protein-Surface Interactions Protein molecules are copolymers made up of defined sequences of the 21 basic L-amino acids. The specific linear sequence of amino acids defines a protein's primary structure. Amino acid side-chains have various sizes, shapes and chemical properties. Some side-chains are acidic and some basic, meaning a protein molecule may vary in charge with changing solution pH. Some amino acid side-chains are relatively hydrophobic and therefore prefer to be positioned so that contact with water is minimized (i.e. within the core of protein). Thus, the primary structure defines a protein's secondary and tertiary structures. Secondary structure involves the formation of or-helices and /2-sheets within a protein, which result from regularly repeating hydrogen-bonding groups in the polypeptide backbone. Association between elements of secondary structure leads to formation of structural domains within the molecule. A three-dimensional tertiary structure then results from association between structural domains, and is stabilized by a complex set of interactions involving hydrophobic forces, electrostatic forces, dispersion forces, and hydrogen bonds [Zubay, 1993]. Tertiary structure is well defined in fibrous proteins (structural proteins) and globular proteins (enzymes, transport proteins). Folding of proteins into these well-ordered structures is entropically unfavourable, and the process occurs spontaneously because of favourable enthalpic interactions. The Gibbs energy of denaturation AN -DG for single domain proteins falls in the range 20-100 kJ/mol, which corresponds roughly to that of 1-8 hydrogen bonds, indicating that the native folded conformation of a protein is only marginally stable [Haynes and Norde, 1994; Privalov, 1979]. Several major driving forces for non-specific protein adsorption have now been identified. Transmission circular dichroism, nuclear magnetic resonance, and microcalorimetry 4 have shown that a significant proportion of a protein's ordered secondary structure can be lost upon adsorption to a surface [Haynes and Norde, 1994]. This suggests that a gain of conformational entropy helps drive the adsorption process. Dehydration of sorbent or protein surfaces has also been shown to influence protein adsorption at a surface [Haynes and Norde, 1994]. Differential scanning microcalorimetry has shown that the total heat capacity Cp of a system often decreases sharply when protein is adsorbed to hydrophobic surfaces. This decrease in heat capacity is consistent with the release of water molecules involved in solvating the surface, a signature of the hydrophobic effect. More hydrophobic proteins therefore show higher maximum adsorption values and less tendency to desorb upon dilution. Electrostatic interactions and redistribution of charge in the interfacial layer also contribute to protein adsorption. Based on these driving forces for protein-surface interactions, neutral, flexible polymer chains attached by one end at the surface should provide resistance to non-specific adsorption. The protective polymer layers are prepared such that long flexible chains extend away from the surface into solution, forming a soft hydrophilic barrier. The hydrophilic nature of the chains can act to reduce any hydrophobic driving force. Also, solute molecules that move into volume that the polymer chains can occupy (excluded volume) decrease the polymer conformational entropy, which is energetically unfavourable. This entropic effect may help to compensate for any gain of protein conformational entropy that adsorption would provide the system. Indeed, it is accepted but not proven in the literature that some combination of these two effects is what leads to reduced protein adsorption on grafted polymer surfaces [Schroen et al, 1995; Gblander etal, 1992; Gombotz etal, 1992; Andrade, 1986]. 5 Indirect studies indicate that entropic repulsion may be the more dominant effect. Schroen et al. [1995] studied lipase adsorption on surfaces covered with adsorbed Pluronic block copolymer (PEOm-PPOn-PEOm; PEO = polyethylene oxide, PPO = polypropylene oxide). When adsorbed to hydrophobic surfaces by the PPO block, the Pluronic formed a highly extended surface layer of PEO tails, and lipase adsorption was reduced to undetectable levels. If the copolymer was adsorbed onto a hydrophilic surface, the PEO tails collapsed onto the surface. Although the amount of copolymer adsorbed still provided a very hydrophilic surface, lipase adsorption was not reduced as significantly. 1.3 Thesis Objective Significant experimental data are now available describing solute interactions with terminally attached polymer chains. Most of this work has focused on trial-based correlation of materials and synthesis methods with protein adsorption or cell adhesion experiments [Merril, 1992; Golander et al, 1992; Gombotz et al., 1992; Desai and Hubbel, 1991, Bergstrom et al., 1992]. Often, however, incomplete characterization of the experimental system leads to ambiguous or incomplete conclusions. It would be very useful to understand these systems well enough to predict a priori how a solute molecule will interact with a particular terminally attached polymer brush. This involves (i) defining how the polymer will change its equilibrium distribution when a solute particle is interacting with it as well as (ii) understanding the energetics accompanying the changes in polymer distribution and changes in contact between different components of the system. The goal of this work is therefore to develop a general model to study the interaction of a solute particle with polymer that is terminally attached at an interface. The first objective is to calculate the polymer density distribution in the presence of a particle whose position may be 6 varied. In this way we hope to gain a better understanding of how the polymer will redistribute itself as the particle penetrates toward the surface. The second objective involves calculation of the interaction energy between the polymer and the particle in order to understand better how the polymer repels a particle. The effect of varying polymer chain length and surface density on the polymer distribution and interaction energy will be investigated. Also the effect of varying particle size and interaction parameters between system components will be investigated. It is hoped that the model will help us to understand how a surface might be protected most effectively from protein adsorption using terminally attached polymer. Also, the effects of system properties on polymer-particle interactions are of interest for other applications such as separations using end-grafted polymer surfaces. 1.4 Terminally Attached Polymers at Interfaces A polymer is a large molecule composed of a number of small, chemically simple subunits. The subunits from which the polymer is formed are called monomers and they define the type of polymer (see Table 1.1). Polymer molecules are often linear sequences of repeat units (e.g. polyethylene), but may also be branched chains and even cross-linked networks of chains forming very large molecules. The degree of polymerization N is defined as the number of repeat units forming a given polymer, and hence the molecular weight of the molecule M is equal to its degree of polymerization multiplied by the repeat unit molecular weight. If all repeat units are identical, a polymer is called a homopolymer; copolymers on the other hand contain a mixture of different monomers. The different types of copolymer include random, alternating, block, and graft copolymers (see Table 1.2). For example, if one monomer type forms a linear chain, while a second monomer type is coupled to an active site along that 7 chain (other than the ends), a graft copolymer is formed. Both block and graft copolymers are often used to modify interfacial properties for reducing protein adsorption. Table 1.1 Several polymers with their monomers and repeat units. Polymer Monomer Repeat Unit polyethylene C H 2 = C H 2 - ( C H 2 C H 2 ) N - polyvinyl chloride CH 2=CHC1 -(CH 2 CHC1) N - polyethylene glycol HOCH2CH2CH2OH - ( C H 2 C H 2 0 ) n - polyvinyl alcohol C H 2 = C H O H - ( C H 2 C H O H ) N - Table 1.2 Schematic of different copolymer structures. (a) Random Copolymer -AAABABBABAABBB- (c) Block Copolymer -AAAABBBBAAAA- (b) Alternating Copolymer (d) Graft Copolymer -ABABABABABABAB- -AAAAAAAAAA- • w w CO 1 When linear polymer molecules are placed in solution they may adopt a great number of conformations whose average is that of a random coil. Although the bond angles between atoms forming the backbone are fixed (e.g., 109.5° for C-C), the energy barrier associated with rotation 8 about those bonds is often on the same order of magnitude as the thermal energy of the solution. As a result, the path from one end of a long polymer molecule to the other changes over time. Thus, describing geometric properties of polymer molecules in solution requires a statistical analysis. For example, the average path and path length of a chain can be determined using a statistical-mechanics formalism known as the random walk [Kuhn, 1934]. In this manner, the average distance between ends of a polymer molecule can be calculated by considering the chain as a sequence of N randomly oriented bonds each of length /. The well known result is [Flory, 1953]: (R2) = l2N (1.1) From this result, the radius of gyration (R2 ^ of a random coil chain is given by (R2g) = ̂ l2N (1.2) which is defined as the root mean squared distance of monomers from the center of mass. It has also been noted that a polymer chain will form a somewhat more expanded coil in a very good solvent, de Gennes (1979, 1980) proved that the coil radius in this case may be given by Rg «IN3'5 (1.3) where numerical constants of order unity are ignored. The central result here is that Rg oc N" where a lies in the range from 1/2 in a theta (poor) solvent to 3/5 in a good solvent. Polymer chains are not completely flexible however, due to the fixed bond angles and steric interactions between side groups. Flory was the first to account for this by multiplying the random walk result (Equation (1.2)) by a stiffness parameter CX, which is determined experimentally. A parameter known as the persistence p may also be used to account for chain stiffness and it is equivalent to Cx/6. Finally, rescaling of N and / can be used to account for chain stiffness [Kuhn, 1934]. For example, if a statistical segment length IK (Kuhn length) is chosen such that lK = bl and the number of segments is chosen as NK = NI c, the average chain dimension becomes (R2) = IRNk = (b21c)l2N. The chain stiffness will be properly accounted for as long as b and c are chosen such that (b21c) = C M =6p [Fleer et al., 1993]. In the following discussion the chain persistence is ignored for simplicity, but can be reintroduced in a straightforward manner without loss of generality. 1.4.1 Distributions of Terminally Attached Chains 1.4.1.1 Experimental Measurements and Scaling Predictions Graft copolymer synthesis is the most common method used to terminally attach polymer at an interface. For example, polystyrene microspheres onto which PEG chains are grafted are commercially available for various applications including oligopeptide synthesis and chromatography [Bayer and Rapp, 1992]. Polymer chains may also be immobilized at a surface by block copolymer adsorption [Van Alstine and Malmsten, 1997; Kent et al., 1995; Schroen et al., 1995]. If one block of a copolymer has a very strong affinity for the surface it will adsorb preferentially leaving the other block freely extended into solution. In this situation, it is possible for adsorbing molecules to laterally displace molecules already at the interface, giving rise to a surface-mobile adsorbed layer (especially when adsorbed at a liquid-liquid interface). A third technique for preparing surfaces with terminally attached chains involves self-assembly of 10 monolayers or bilayers using polymer-derivatized lipids [Du et al., 1997; Woodle and Lasic, 1992]. Restriction of one end of a polymer molecule to an interface leads to very distinct chain configurations and energetics. The degree of surface coverage a (chains/unit area), or distance d between grafting points (d & a ~U2\ in large part determines the grafted chain conformation. A crude picture of chain conformation can therefore be obtained by calculating the reduced surface coverage a* = n (R2^JY,, which is approximately the ratio of an undisturbed chain's cross- sectional area to the average surface area per grafted chain, where S=l/cr (area per chain) [Baranowski and Whitmore, 1995]. If the grafting points are farther apart than the chain radius of gyration, a random coil configuration, commonly termed a "mushroom", is adopted by the chains. In the mushroom regime the reduced surface coverage is less than unity and the thickness of the polymer layer h is independent of surface coverage and will increase with N V , where v= 0.5 to 0.6 in a poor to good solvent respectively (corresponding to Equations (1.2) and (1.3)). Within the mushroom regime, grafted-chain conformations are extremely sensitive to enthalpic interactions, commonly expressed in terms of a binary Flory interaction parameter x<y For example, if a favourable energetic interaction exists between chain segments and the surface, coil collapse can occur forming a pancake structure on the surface (see Figure 1.1 (d)) [Fleer et al., 1993]. In regimes of surface coverage where cr* >1, chain overlap occurs and the chains extend away from the surface due to excluded volume interactions with one another. A polymer layer under these conditions is often termed a "brush" due to the physical picture of partially ordered, stretched, linear molecules extending from a surface. This regime was first described by 11 Alexander [1977] and de Gennes [1976], who predicted that the polymer configuration minimizing entropically unfavourable chain stretching and excluded volume interactions between chains would result in an equilibrium brush with height h. The equilibrium brush height was found by minimizing an approximate expression for the brush energy derived either using polymer scaling theory [Alexander, 1977] or global energy balance arguments [Fleer et al., 1993]: AA&kT 3h 2 2N17 + vN (1.4) Here, the first term in square brackets represents the energetic cost to stretch a polymer from a Gaussian coil having entropy of 3k/2 per segment or AGk/2 per chain, where k is Boltzmann's constant. The stretching energy is proportional to distance squared (i.e. Hookian), and the ratio h21 (NI)2 gives the fraction of total chain conformational entropy lost due to stretching, where all entropy is lost at the fully stretched length of NI [Fleer et al., 1993]. The second term represents the excluded volume interaction energy between polymer segments, where v represents an excluded volume parameter (equal to \-2%n : see Section 2.1) and Na I31h is equivalent to the average volume fraction <D of polymer in the brush. Ignoring numerical constants of order unity, the equilibrium (AA = 0 ) height is therefore given by [Milner, 1991] h^Nl(val2)U3 (1.5) which indicates that the brush height scales linearly with N. As suggested by Equation (1.2) the unstretched chain dimension Rg therefore scales with Nm. Equation (1.5) also suggests that brush height increases with <r1/3, so that the maximum volume fraction increases with a2'3. 12 However, two major assumptions of this model are that (i) all free chain ends lie on the plane at the outer edge of the brush and (ii) the volume fraction of polymer in the brush assumes a constant density step profile (Figure 1.1 (a)). Both assumptions are at best approximate. Thus, the energy expression overestimates both the stretching cost and the excluded volume energy, but cancellation of errors results in the physics of the system being described reasonably well [Milner, 1991]. Following these early model results, experimental brush density distributions have been measured using small angle neutron scattering (SANS) and neutron reflectivity to provide a much clearer picture of chain conformations. Cosgrove et al. [1987b] measured the density profiles of PEO grafted on nonadsorbing, deuterated polystyrene using SANS. The polymer density profile exhibited a clear maximum several nanometers from the surface and then decreased slowly (parabolically) away from the surface. At the outer edge of the brush, the polymer volume fraction was seen to decrease very slowly giving rise to a long tail region (Figure 1.1 (c)). Density profiles for PS-PEO adsorbed from d-toluene onto quartz measured by Field et al. [1992] using neutron reflectometry were also of the form shown in Figure 1.1 (c). Similar profiles were measured by Cosgrove et al. [1991] using neutron reflectivity, for polystyrene-poly(2-vinypyridine) (PS-P2VP) block copolymer adsorbed onto mica from toluene. In this solvent, the PS forms a highly extended layer, while the P2VP remains strongly adsorbed. Auroy et al. [1991a, 1991b] studied the structure of end-grafted polydimethyl siloxane (PDMS) on silica using SANS and found that in a poor solvent the polymer layer was collapsed and gave good agreement with the step-profile model. However, in a good solvent the layer was highly extended such that such that the volume fraction Q>(z) decreased very slowly at high z. Auroy et 13 al. [1991a] did not see a depletion layer next to the grafting surface but they suggest that there may be an attractive interaction between the PDMS polymer and the surface (Figure 1.1 (d)). Z Z Figure 1.1 Several possible polymer density profiles are shown for terminally attached polymer chains, (a) The Alexander-de Gennes brush assumes a step profile, (b) Analytical SCF theory predicts a parabolic brush, (c) Real terminally attached polymer will exhibit a maximum slightly away from the surface, which advanced numerical models also predict, (d) A strongly adsorbing surface can lead to a pancake structure. Cosgrove and Ryan [1990b] successfully measured the density profile of end-adsorbed methacrylate-terminated polyethylene glycol methyl ethers (Mn=5000 g/mol) on polystyrene microparticles using SANS. They found that the density of polymer was a maximum at the surface and decayed rapidly away from the surface. This was attributed to low surface coverages (1-4 mg/m2), and the fact that the surface was strongly charged and an adsorbing surface for the PEG chains. Cosgrove [1990a] also presents data for PS grafted to silica in which there are 14 different polymer-surface interactions. Here the polymer-surface interaction was characterized by a dimensionless adsorption energy %s = {u° -u^/kT that represents the difference between the energy of a solvent-surface contact u" and a polymer segment-surface contact u2 [Silberberg, 1968]. For high %s the segment density was found to be a rapidly decaying function away from the surface (Figure 1.1 (d)); for low %s a parabolic-type profile was seen with a depletion region next to the surface (Figure 1.1 (c)). At very high brush densities, however, segment adsorption can result in no depletion layer being observed [Auroy et al, 1991a]. Kent et al. [1992; 1995] carried out experiments in which the chain length and surface density of a brush were varied independently over a large range. Langmuir monolayers of polydimethylsiloxane-polystyrene (PDMS-PS) copolymer were prepared at an ethyl benzoate-air interface, and density profiles were measured using neutron reflectivity. By varying the surface pressure n it was possible to vary the density of chains in the brush and thereby test the accuracy of various analytical models. In contrast to the scaling predictions, Kent et al. found that brush height scaled as hxa°'22N0M. These results suggest that the analytical models are often limited by their inherent assumptions and are not truly applicable to the majority of experimental brushes. The distribution of terminally attached polymer chains in a poor solvent is more poorly understood because the system cannot undergo ordinary phase separation. Uchida and Ikada [1997] imaged the surface of poly(2-dimethylamino)ethyl methacrylate (PDMAE-MA) brushes grafted on polyethylene terephthalate films in water using scanning atomic force microscopy. They observed surface inhomogeneities, which were measured reproducibly and attributed to the fact that the polymer was below its lower critical solution temperature in water. They also found that segment clusters would tilt in the direction of the AFM scan due to friction between the 15 AFM tip and the chains, with less load on the tip being required to induce tilting in longer chains. This suggests that clumping of polymer segments can occur in a very poor solvent along with brush collapse. 1.4.1.2 A d v a n c e d Brush Models Milner, Witten and Cates [1988a, 1988b] and Zhulina, Pryamitsyn and Borisov [1989] independently developed analytical self-consistent-field models to describe brush distributions that provide essentially identical results. The basis of their model is the assumption that the set of all possible chain conformations may be replaced by the average chain trajectory, neglecting fluctuations about this average, when the chains are strongly stretched as in a brush [Semenov, 1985]. The main result is that the brush density profile is found to be a gradually decreasing function in the z-direction (normal to the surface). It is expressed in analytical form as a parabola as shown in Figure 1.1 (b), which is given by 0>(z) = (B/v)[h2 -z2] (1.6) where h is the equilibrium height given by h-(\2va/7r2y/3N and B is a prefactor equal to TT2/SN2 [Milner, 1988a]. The brush density profile is predicted by Equation (1.6) to be a maximum at the surface and also to terminate abruptly at distance h from the surface. Similar to the Alexander-de Gennes brush, the equilibrium height is found to scale with N and c 1 / 3 . In this case, however, the brush is slightly more extended than the step-profile brush (h/hstep=l.3). Moreover the analytical SCF theory predicts a finite probability to find free end segments throughout the brush. Although many physical systems of interest may not lie in the regime of strong stretching, the analytical result is exact within this assumption and gives simple expressions to describe brush properties. Because of the compact analytical form of the step- 16 profile brush and parabolic brush models, they have received a great deal of attention [Kenworthy et al, 1995; Baranowski and Whitmore, 1995; Milner, 1991]. More recently, advanced simulation techniques and lattice theories have been applied to the determination of brush distributions. Detailed computer experiments including numerical SCF calculations and Monte Carlo and molecular dynamics simulations have all shown good agreement with measured density profiles [Baranowski and Whitmore, 1995; Lai and Binder, 1991; Chakrabarti and Toral, 1990; Murat and Grest, 1989; Cosgrove et al, 1987a]. Both approaches are able to predict a depletion layer next to the surface leading to a maximum in the density, followed by a region of decreasing density that terminates in a smooth tail region. The brush extension has been found to depend on the solvent conditions and surface adsorption. Computer experiments have also been useful for exploring system conditions that are not easily obtained in experiments, such as dependence of brush conformation on chain length and surface density. In addition to their use as a test of model accuracy, simulations (both Monte Carlo and molecular dynamics), have also provided an accurate, detailed picture of tethered chain conformations and their influence on segment density distributions. Monte Carlo simulations, for instance, have been used by Chakrabarti and Toral [1990] and Lai and Binder [1991] to determine the true scaling behaviour of brush height and maximum segment density. Molecular dynamics calculations addressing the same questions have also been performed by Murat and Grest [1989]. Simulated chain lengths range from N = 20 to 99 and surface fractions from s = 0.04 to 0.2, thereby providing an accurate picture of brush structure over the realistic range of brush synthesis conditions. 17 Wijmans et al. [1992] carried out a detailed comparison of brush distributions calculated using an advanced numerical SCF model (Chapter 2), which makes no a priori assumptions about chain or brush distributions, with those calculated by the analytical theory of Zhulina et al. [1989]. The analytical theory cannot predict a depletion zone next to the surface; it also predicts a sharp end-point of the brush due to the central model result of a parabolic segment-density profile. As a result, important features of real brushes, which are accurately captured by numerical SCF theory, are lost in the analytical model. Nevertheless, the two analytical models have received considerable attention due to their simplicity. Moreover, both analytical models predict the exact behaviour of the brush in the limit N —> oo or under fully stretched chain conditions, which should be approached, at least approximately, in an athermal (i.e., a good) solvent. Lai and Binder [1992] compared scaling prediction from analytical models to results from Monte Carlo simulations for low-molecular-weight brush density profiles in near theta solvent conditions over the surface fraction range of o = 0.075 to 0.125. Simulated density profiles were characterized by a pronounced depletion zone near the grafting surface and a smooth tail region, similar to that shown in Figure 1.1 (c). Thus, even under good solvent conditions, the analytical models fail to predict important brush characteristics. Nevertheless, normalization (0(z)/cr1/2 vs. z I Na1'2) of simulated density profiles yields a universal curve, which is accurately predicted by scaling results except near the brush boundaries. Under poor solvent conditions, the accuracy of the analytical models is significantly worse. Grest and Murat [1993] carried out molecular dynamics simulations of brushes in poor solvent and found that the parabolic density profile predicted by analytical SCF theory was in significant error. Unlike free chains, the tethered chains in a brush cannot phase separate under 18 poor solvent conditions. Using Monte Carlo simulations, Lai and binder [1992] were the first to observe lateral phase inhomogeneity of brushes in a poor solvent. The simulated chains formed collapsed clusters of segments when the solvent quality became poor enough to significantly favour polymer-polymer contacts. Grest and Murat [1993] also saw this lateral phase separation in their molecular dynamics simulations. At low enough surface coverage (a = 0.03), 50 and 100 segment chains were seen to collapse into clusters leaving bare patches or holes through which the grafting surface was exposed. They saw a decrease in bare surface with increases in chain length and predicted that the holes would disappear with long enough chains, perhaps forming a dimpled brush surface. Yeung et al. [1993] studied the effect of solvent quality on brush phase stability using a numerical mean-field model combined with random phase approximation. They also concluded that as solvent quality decreased, the brush became unstable and would form a dimpled surface, in agreement with the AFM work of Uchida and Ikada[1997]. 1.4.1.3 B r u s h e s with Free Polymer in Solution By interacting with the solvent, grafted chains and surface, free polymer in solution can influence brush configuration. Lee and Kent [1997] used neutron reflectivity to measure the density profiles of PDMS-PS copolymers adsorbed at the interface between air and ethyl benzoate solution containing free polystyrene. Selective deuteration of the adsorbed polystyrene (N - 1625 repeat units) and free polystyrene (N = 413 and 3846) in solution allowed them to measure the density profiles of each component independently. They observed a more compressed brush density distribution in the presence of the free polymer, and the tail region became less pronounced with higher molecular weight free polymer. They also found that there was significant penetration of the free polymer into the brush, even at a * =12, which increased 19 with decreasing brush surface density. The degree of penetration was also seen to increase significantly with decreasing free chain length. Comparison of their experimental results to predictions of a modified analytical SCF theory showed that the free polymer penetrated much more than the model predicted [Lee and Kent, 1997]. Numerical SCF models have also been used to study the structure of the brush in the presence of free polymer. Wijmans and Factor [1996] carried out calculations, which simulated the experimental conditions of Lee et al. [1994] and found good agreement in terms of both density distributions and brush height. They were also able to reproduce the trends of brash height and penetration with varied surface density and free polymer chain length. They suggest that the decrease in brush height is a result of a shift in the balance between chain stretching and excluded volume interactions between brush chains. Penetration of the free polymer acts to screen excluded volume interactions in the brash, thereby reducing its height. However, the model consistently underestimated the brash height in the presence of free polymer. The model also underpredicts the penetration of the free polymer into the brash, although not as badly as analytical SCF models. 1.4.1.4 Brush Distributions - Summary From this earlier work, a picture of polymer brashes in solution has emerged, in which long chain molecules stretch away from an interface, forming a partially ordered layer. This results from a balance between the entropic penalty associated with chain stretching and the unfavourable excluded volume interactions between neighbouring chains. Often, a depletion zone is observed near the surface, due to the entropic penalty imposed on chains by the impenetrable grafting surface. Unfavourable interaction between the polymer and solvent may lead to a more compressed layer, since the excluded volume term is less significant. In lower 20 ranges of surface coverage where the chains form a brush in a good solvent, decreased solvent quality may lead to lateral phase separation and thus exposure of the underlying surface. Also, a favourable segment adsorption energy may lead to absence of the depletion zone in a dense brush, or a rapidly decaying density profile in a less dense brush (pancake). Free polymer in solution will compress the brush, due to the balance between osmotic pressure of the brush and bulk solution, and this compression increases with increasing free polymer concentration. Also, significant free chain penetration into the brush is observed, which is greater for shorter free chains, but decreases with increased brush density. Numerical and analytical models are also very useful for understanding (visualizing) and predicting behaviour of polymer brushes. Overall, the analytical models have shown value because of their simple form and ability to capture the basic physics of brushes. They do however contain some unrealistic assumptions, such as the step-profile brush model, which yield an unrealistic picture except in the case of a very high density brush. Also the strong stretching assumption for the parabolic brush model has been shown poor in most regimes of brush density and chain length that are experimentally obtainable. More complex computer models, such as numerical SCF, Monte Carlo and molecular dynamics models provide a significantly better picture of brush density profiles. Molecular dynamics simulations give presumably the best representation of a real system. However, due to the large computer capacity necessary for both MD and MC simulations, numerical SCF models have been widely used, showing good agreement with results of experiments and the more complex computer simulations. Typical computing time for SCF calculation of brush distributions is on the order of minutes for a desktop PC, which makes them convenient for investigation of a wide range of brush parameters. 21 1.4.2 Compression of Terminally Attached Polymers Adsorbed polymer layers are known to prevent flocculation of colloidal particles in solution [Napper, 1983]. Early attempts to understand this phenomenon and identify conditions that would result in colloid stabilization were carried out by Hesselink [1971], Meier [1967], and Clayfield and Lumb [1966]. They used Monte Carlo simulations and models to calculate interaction potentials between particles on which polymeric stabilizers were adsorbed. Experimental verification of these predictions was provided by Dorozkowski and Lambourne [1971], who measured the repulsive force between colloidal particles experimentally using a surface balance. The repulsion between adsorbed layers of polymer was believed to result from loss of chain configurations and also a change in the free energy of mixing of polymer and solvent [Meier, 1967]. More advanced experimental equipment has allowed the force between two brushes to be measured directly. This has been carried out using the surface force apparatus (SFA) and the atomic force microscope (AFM). Reviews by Luckham [1996] and Claesson et al. [1996] discuss recent advances in the technology. The SFA employs two curved mica surfaces (convex) placed such that their cylindrical axes are perpendicular to one another. Polymer may be adsorbed or terminally attached to the mica surfaces, which are then moved together in small increments (nm scale) using a piezoelectric transducer. The force between the mica sheets is measured as a function of separation distance. The AFM uses a small cantilever with a very fine tip, which contacts a surface at positions controlled by a piezoelectric transducer. In very simple terms, the bending of the cantilever, as the surface is raised or lowered, is measured using a laser reflected off the back of the cantilever into a split photodiode. By knowing the spring constant 22 of the cantilever it is possible to calculate the applied force as a function of distance from the surface. Hadziioannou et al. [1986] measured the force between surfaces of adsorbed poly(2-vinyl pyridine)-polystyrene (P2VP-PS) block copolymers using the SFA. They studied the case of PS in a good solvent (toluene) and found that the onset of repulsion occurred at about ten times the radius of gyration of a free PS chain. This was in contrast to adsorbed homopolymer layers where a repulsion onset was seen to occur at about 3RS, suggesting that a highly extended polymer layer was formed by the block copolymer under good solvent conditions. They also found that the force between surfaces increased monotonically with decreased separation, with the onset of repulsion occurring at higher separations for longer PS chains. Similar results showing a monotonically increasing repulsive force with compression have been shown by others [Kenworthy et al., 1995; Taunton et al., 1990; Taunton et al., 1988; Patel et al., 1988]. Hadziioannou et al. [1986] also observed a weak attraction between brushes in a less- than-theta solvent (cyclohexane at 26 C). This was attributed to net favourable polymer-polymer contacts in the poor solvent. Upon higher compression, the repulsive force increased more rapidly than in a good solvent, due to the more collapsed nature of the brush. Numerical SCF calculations by Van Lent et al. [1990] suggest that an attraction between brushes will be felt if the solvent does not meet theoretical theta conditions (xu = 0.5). They attribute this attraction to the fact that terminally attached chains have essentially zero translational entropy and thus their phase behaviour is similar to that of infinite molecular weight chains. Various attempts have been made to correlate experimental brush height data with scaling behaviour. This is difficult for two reasons. First, surface density and chain length are difficult to vary independently during brush synthesis. As a result, a clear dependence of brush height on 23 N at constant o does not emerge from experiment, making theoretical comparisons difficult. Second, brush height cannot be measured directly, but instead is generally taken as the separation distance where the onset of brush repulsion occurs. It is unclear how this measured property relates to the actual brush height, and if the two properties scale linearly. This latter point has generated much confusion in the literature. Hadziioannou et al. [1986], for instance, report that brush height h scales linearly with TV, while Ansarifar and Luckham [1988] and Taunton et al. [1990] reported a dependence of JV07 and N06, respectively. Scaling theory [de Gennes, 1985], which assumes a step-profile in density, predicts that the force to compress two layers of terminally-attached polymer whose surface coverage lies in the mushroom regime is given by In Equations (1.7) and (1.8) numerical prefactors are ignored, and D represents the distance between the two grafting surfaces. In Equation (1.8) the first term in brackets gives the repulsive force due to the increase in osmotic pressure in the brush (higher excluded volume), while the second term accounts for the energetically favourable decrease in chain stretching. In the limit of close compression, where D«h, the osmotic term dominates and the force is predicted to scale with D'9/4. (1.7) and that the compression force between two brushes on planar surfaces is given by: (1.8) 24 The energy to compress a brush with a parabolic segment density profile has also been determined and is given by [Milner, 1988a]: A(u) = N(a vy'1 (n2 IT2) v l / 3 1 1 - — + -W 2u 2 10 -w (1.9) Here u is given by h'/h, where h' is the new brush height due to compression. Due to the "softer" parabolic density profile, the brush is predicted to have a lower compression energy at small compressions. Both models assume complete compression of the terminally attached chains, which excludes (experimentally relevant) cases where there is interpenetration of brushes. Ansifar and Luckham [1988] compared scaling-theory predictions to their experimental results for the compression force of poly(2-vinylpyridine)-poly(t-butylstyrene) (P2VP-PBS) brushes adsorbed in toluene. They found that scaling theory underestimated the rate of increase in force with compression, which was found to be proportional to D~4. In contrast, Taunton et al. [1990, 1988] measured the force between two brushes of PS in toluene under compression and found that the de Gennes scaling dependence of force provided a good representation of the data. Since the scaling results ignore numerical prefactors, it was necessary to adjust the curve with a single constant whose value was found to be about 1.5. Taunton et al. [1990] also compared the measured compression force with predictions of the model of Milner et al. [1988a]. The force calculated from the model was in very good agreement with experiment if the surface density was set to a value slightly higher than the experimental value. Patel et al. [1988] have derived a more general expression for the brush compression energy, based on the de Gennes result, which introduces a parameter describing the solvent quality and numerical prefactors {k\ and kq) describing the contribution of the osmotic and 25 stretching terms to the energy, respectively. Although their expression is more complex, similar results are obtained. A detailed study by Kenworthy et al. [1995] investigated the pressure-distance relationships of lipid bilayers containing PEG-derivatized lipids. The distance between liposome bilayers as a function of osmotic pressure was measured using X-ray diffraction. Bilayers with different PEG chain lengths and different mole fractions of derivatized lipid were constructed and used to study the effect of varied chain length and surface density on brush-brush interactions. They characterized the different regimes of surface coverage as (i) interdigitating mushrooms (where mushrooms could fit into gaps between mushrooms on the opposite surface), (ii) mushrooms and (iii) brushes. For all cases, they observed a monotonic increase in osmotic pressure with decreased separation. They also observed an earlier onset of repulsion and higher repulsive pressures with increased N or <j. By comparing various models to their result, Kenworthy et al. [1995] found that compression of mushrooms in the interdigitating regime was best modeled assuming compression with no interdigitation, since pressure was underestimated by all models at large separations. They reached this conclusion even though at higher compressions it appeared that some interdigitation did occur. In the brush regime, Kenworthy et al. found that the de Gennes and Milner results gave good agreement at strong compression where the osmotic term dominates, but both models underestimated the repulsive pressure at weak compressions. They gave no explanation for the discrepancy, but the most likely explanation is that the analytical models underestimated the brush height, due to the sharp termination of the step or parabolic density profiles. 26 The analytical models of de Gennes and Milner also suffer from the assumption that no interpenetration occurs between the opposing layers as they move toward one another. Moreover, all molecular detail is lost in these models since they are based on energy balance arguments. Both of these limitations can be removed through more powerful computer simulation techniques. Murat and Grest [1991], for example, investigated the interaction of two brushes using molecular dynamics simulations and found that a significant degree of interpenetration occurs. They were able to show good agreement between their calculated force- distance curves and experimental results. Chakrabarti et al. [1994] confirmed that significant interpenetration occurs using Monte Carlo simulations. Brush-brush interactions have also been investigated using numerical SCF models. Wijmans et al. [1994a] used numerical SCF theory to accurately calculate force-distance curves at weak compression and showed that the analytical model of Zhulina et al. [1989] underestimated the force at weak compression due to the absence of the tail region in the parabolic density profile. Martin and Wang [1995], also using numerical SCF, confirmed the importance of the tail region for the onset of repulsion. A more advanced numerical SCF model, which takes bond neighbour correlations into account, has recently been developed by Ruckenstein and Li [1997]. Calculated force-distance curves are in good agreement with experiment, in contrast to analytical models, which show large deviations in the regime of low compression. The numerical model required no adjustable parameters, although one could adjust the polymer-solvent interaction parameter. Ruckenstein and Li [1997] investigated the importance of chain interpenetration between the brushes and showed that the compression force was overestimated by up to an order of magnitude when interpenetration was not allowed. 27 1.4.3 Particle Interactions with Terminally Attached Polymers The interaction of finite-sized particles with terminally attached polymers has received considerably less attention than the interaction between two end-grafted polymer surfaces, at least on the fundamental level. This is due to the difficulty of making experimentally meaningful measurements and the difficulty in modeling these systems. The majority of experimental work in this area involves protein adsorption experiments on surfaces modified with terminally attached polymers [Du et al, 1997; Gombotz, 1992; Gblander, 1992]. Some useful information is also now being provided by chromatographic separation of proteins using terminally attached polymer [Brooks and Miiller, 1996]. 1.4.3.1 Experimental S t u d i e s Experimental work describing the interaction of proteins with end-grafted polymer does not generally provide solid conclusions about how polymer properties influence the interactions. Different optimum chain lengths of end-grafted PEG have been reported for minimizing protein adsorption on PEG-modified surfaces. For PEG hydrogels, Gblander et al [1992] reported that adsorption decreased with increasing polymer molecular weights up to 1.5 kDa, above which only a marginal decrease is observed. In attempting to minimize adsorption of fibrinogen and bovine serum albumin on polyethylene terephthalate films, Gombotz et al [1991] reported an optimum PEG graft length of 3.5 kDa, above which little decrease in adsorption was seen. Desai and Hubbel [1991] reported that adsorption of albumin and fibrinogen onto PEG-grafted polyethylene terephthalate surfaces was minimized when PEG chains of molecular weight 18 500 g/mol or higher were used. 28 The difficulty in drawing solid conclusions about the optimum properties of grafted chains is due in large part to the interplay between surface density and degree of polymerization. Characterization of surfaces using electron spectroscopy for chemical analysis (ESCA) [Gombotz et al., 1991; Desai and Hubbel, 1991] showed that higher molecular weight PEG was present on the surfaces in lower quantities. It is believed that the experimentally achieved surface densities are often lower for large polymer chains since they have a higher excluded volume and lower diffusion to the interface during synthesis. Thus, a given chain length may appear to give optimum protection because it resulted by chance in the best combination of chain length and surface density for protein exclusion after experimental synthesis. Du et al. [1997] studied the behaviour of PEG derivatized lipid bilayers with respect to protein adsorption and cell adhesion. Like the work of Kenworthy et al. [1995], they were able to vary the PEG chain length (M = 750 - 5000) and surface density (0-5 mol% PEG-lipid) independently. For the range of PEG-lipid substitution used, the longest chains were seen to give the best reduction in protein and cell adhesion. They also observed a lower degree of protein adsorption and cell adhesion with increasing surface density of the PEG chains. For M= 5000 chains, the most significant decrease in adsorption and adhesion was seen up to about 1 mol% PEG-lipid, which is just above what they calculated as the theoretical point where complete surface coverage is achieved (0.7 mol%). Due to the fact that some protein still adsorbed after theoretically complete surface coverage, however, they concluded that complete surface coverage was never truly obtained. It appeared that the PEG chains were not randomly distributed, but instead arranged in denser and less dense regions of coverage. Yang and Yu [1997] prepared self-assembled monolayers of end-functionalized PEG chains on glass slides and studied the adsorption and translational diffusion of bovine serum 29 albumin (BSA) on the surfaces. Using a fluorescent dye, they found that the level of adsorption decreased to about 6% of that on the bare glass surface. They also found that the surface diffusion coefficient measured using fluorescence recovery after photobleaching (FRAP) increased by almost an order of magnitude on the coated surfaces. Complete recovery of fluorescence was not observed, suggesting that some of the BSA was not surface-mobile. In contrast to adsorption on bare silica, however, adsorption was completely reversibly. Yang and Yu [1997] also imaged protein sorbates on PEG covered and bare surfaces using an AFM. They observed regular ellipsoid patterns on the PEG-covered surface, which had the correct dimensions of BSA molecules, while the bare surfaces showed no regular patterns. From this information, they speculated that PEG reduced the magnitude of interactions between the protein and surface to a level where native-state conformations were retained. 1.4.3.2 Modeling of Brush-Particle Interactions The empirical model of Joen et al. [1991a; 1991b] represents the first attempt to analyze the interaction of a spherical solute with a brush, in this case end-grafted polyethylene oxide (PEO). Interaction energies were calculated as the sum of a steric repulsion term, a van der Waals attraction term and a hydrophobic interaction term. The steric repulsion term is based on the brush model of de Gennes and makes use of the compression energy expression derived by Patel et al. [1988]. Thus, the model assumes the brush is completely compressed with no splaying as the particle approaches. It predicts that solute repulsion increases with increasing particle penetration and with increasing chain length at a given chain density. It also predicts, somewhat surprisingly, that the optimum chain density for solute repulsion depends upon the particle size. For larger proteins of radius R = 60-80A, the optimum distance between grafting 30 points d was found to be slightly more than that for small proteins of radius R = 20A (d = 13- 17Avs. 9-11 A). A molecular dynamics simulation of interactions between lysozyme and end-grafted PEG chains by Lim and Herron [1992] showed that the mobility of the chains may be important for reducing protein adsorption. Their simulations assumed a weak attraction between the polymer and some surface regions of the protein. In the simulations, the protein always diffused away from the surface although weak attachment of polymer segments to the protein sometimes occurred. The polymer chains were able to move with the protein until the contacts were broken. Lim and Herron [1992] also calculated the energy required for the protein to move into the surface of the brush and found a slight attraction at very weak penetration (4A) followed by strong repulsion. In some cases there was a maximum in repulsion energy at very weak penetration that they attribute to the protein squeezing through a small opening in the PEG chains. In order to treat brush interactions with a small particle, Subramanian et al. [1995] have modeled the interaction of a finite-sized cylinder with a brush using the equation for brush compression of Milner et al. [1988a]. In this way, they consider the case where only brush compression occurs and show that the repulsive force is always increasing with approach of the particle. They ignore the case where polymer chains may spread around the particle, although they acknowledge that this would be important for small particles. One clue as to the importance of particle size is offered by the numerical SCF work of Wijmans et al. [1994b] who studied the interaction potential of two spherical particles of radius ranging from 2 to 10 lattice units (on the size scale of proteins) whose surfaces were covered with grafted chains. They found that the repulsion between particles was much lower than that predicted using the Derjaguin 31 approximation for repulsion between two flat surfaces. This was attributed to the ability of polymer chains (N= 50) to move out of the gap between the approaching particles. This kind of molecular detail is required to see the correct interactions between polymer and very small particles. Recently, Murat and Grest [1996] reported MD simulation results for an AFM tip interacting with a brush, from which some insights for an interacting particle can be inferred. The simulations show that the applied force required increases with brush compression. The increase in force is less than that predicted by the simple model of Subramanian et al. because the chain segments are in fact able to escape from under the AFM tip. Larger tip sizes lead to a more rapid increase in compression force since it is more difficult for the chains to escape. However, since the chain segments are able to rearrange themselves around the tip, the segment density between the tip and the surface was found to remain approximately constant (the same as that of the undisturbed brush), in sharp contrast to predictions of analytical models. 32 2 Model Outline 2.1 Origins of Nonideal Mixing in Liquid Solutions In order to understand the properties of a solution, it is often useful to compare behaviour of a mixture to that of an ideal solution as defined by Raoult's law. The activity ai of any component /' is then equal to its mole fraction xt. Deviations from Raoult's law are observed when this condition is not met and mixing therefore results in a nonzero change in enthalpy (AHm * 0) or nonzero excess entropy (SB ^0). Under such conditions, the ratio aj Ixt is no longer equal to unity, but is instead set equal to y\, the activity coefficient of /'. Mixing nonidealities are known to result from asymmetries in the energetics of the molecular interactions occurring in the mixture. For instance, water must break energetically favourable hydrogen bonds with itself in order to solvate a non-hydrogen-bonding solute. As a result, AHm is nonzero and y\ is often far from unity for the dilute solute component. In polymer solutions, mixing leads to large negative values of SE due to the connected nature of the polymer chain, which severely limits the number of ways in which the polymer can be arranged in the mixture. As a result, nonideality is observed even though AHm is near zero. Moreover, nonidealities occur at extremely small mole fractions of polymer. This latter effect arises because of the large size difference between polymer and solvent. Thus, the mole fraction of polymer remains small even as its weight fraction approaches 0.5. The unusual mixing properties of polymer solutions were first recognized in a formal way by Flory [1953; 1942] and Huggins [1942], whose combined theory has served as a basis for understanding polymer-mixing thermodynamics for over 50 years. 33 2.1.1 Entropy of Mixing Polymer in Solvent In order to conform to ideal solution behaviour, a solution must meet two requirements [Flory, 1953]. The first is that the entropy of mixing is equal to = - / £ > , In x, (2.1) i where R (J/K-mol) is the gas constant and x, is the mole fraction of component i. The second is that all mixing is athermal or the heat of mixing AHm is zero. Since the mole fraction of each component must be less than unity, Equation (2.1) is always positive. The increase in entropy is due to the increased number of possible arrangements of the molecules in the system. Modeling a binary mixture using a simple lattice allows the number of different possible arrangements of the molecules to be counted, at least approximately. Flory [1953] showed for a lattice of n0 sites, which contain solvent and n2 solute molecules, the number of ways of arranging them is n l Q = - ^ - (2.2) nx\n2\ which is simply the number of ways of arranging n0 things nx at a time. Through Boltzmann theory AS =-A: Inn (2.3) which reduces to Equation (2.1) after substitution of Equation (2.2). For the case of a polymer solution, Equation (2.1) greatly overestimates the entropy of mixing. This is due to the fact that n2 polymer molecules of degree of polymerization N will have less possible mixing arrangements with n, solvent molecules than an equal volume of 34 solvent-sized solute molecules {i.e. Nn2) would. Flory-Huggins theory treats the polymer molecules as a series of connected segments, and calculates the number of possible ways of inserting the adjacent (connected) segments into the lattice. For example, the number of sites available for placement of molecule i+l is given by v , + 1 =(„ 0 -M)z (z - i r ( l - / , r (2.4) where /' is the number of previously placed molecules, z is the lattice coordination number (i.e. the number of nearest neighbours), and fi is the probability that a site adjacent to the previously vacant one (where the previous segment was placed) is occupied. This probability of occupancy ft in the mean field approximation is simply given by 7V7 / n0. The total number of possible ways of arranging n2 sets of N consecutively adjacent lattice sites correcting for the fact that the molecules are identical is: 1 "2 n0\ , \«2(W-1) 2-1 V " o J (2.5) Substituting this expression into Equation (2.3) and subtracting the reference state of unmixed, amorphous components gives the surprisingly simple expression A S m = - « 0 J R £ ^ l n $ , (2.6) for a mixture of / different components. For very large molecules, Equation (2.6) states the entropy of mixing no longer depends directly upon the mole fraction, but rather on the volume fraction <D . (= Nini I n0) of each component. Figure 2.1 shows a comparison of the entropy of mixing calculated using Equations (2.1) and (2.6). The Flory-Huggins A5"m is approximately 35 half AS'm . It should be noted, however, that Flory Huggins theory slightly overestimate the decrease in mixing entropy of polymer solutions. Thus the true entropy of mixing lies somewhere between the two curves shown in Figure 2.1 [Johansson et al., 1997]. x. or <D. Figure 2.1: Comparison between ideal entropy of mixing and combinatorial entropy of mixing given by Flory-Huggins theory for a polymer chain of 1000 segments. 2.1.2 Enthalpy and Gibbs Energy of Mixing In order to account for enthalpic interactions between components in solution, Flory assumed that only nearest neighbour contacts are significant. Thus, interactions between components that are not immediate neighbours are viewed as being sufficiently weak to make a negligible contribution to the total mixing energy. As with the entropy, a change in energy with 36 mixing is calculated using the unmixed components as a reference state. This is carried out by considering the difference in energy between the like contacts, [1,1] and [2,2], and the unlike contact [1,2] (for the simple case of a binary mixture). Let con, co22 and au be the energy values associated with the respective contacts. It is then possible to write an interchange energy, which is defined by statistical mechanics, as Acou - con (2.7) based on the stoichiometry of ̂ [l, 1] + V2[2,2] = [1,2]. The heat of mixing is then given by AHm =zA<y12n1<D2 (2.8) where component 1 is assumed to be a monomeric solvent. This result is often further simplified by defining a dimensionless parameter Xw which incorporates the interchange energy and changes its basis from the energy change with respect to one contact, to that for interchange of one lattice site. This parameter, called the Flory-Huggins interaction parameter, is defined as %n=zA(on/kT (2.9) and can be measured experimentally; typical values lie in the range of -1 to 1 kT [Polymer Handbook]. Flory's expressions for the entropy and enthalpy of mixing provide a simple, general expression for the Gibbs energy of mixing polymer and solvent cp ™ (2.10) 37 where m is the total number of components. Equation (2.10) can be used to calculate any thermodynamic property of a polymer solution, such as the chemical potential of the solvent or polymer. 2.2 Self-Consistent Mean-Field Theory of Scheutjens and Fleer Flory-Huggins theory deals with an entire solution volume using the mean-field approximation. It is possible to get more detailed information from a lattice model by defining smaller regions in which properties are averaged. In order to study the adsorption of polymer at an interface, Scheutjens and Fleer [1979;1980] treated each layer of a Cartesian lattice with mean-field averaging as shown in Figure 2.2 (a). They developed a self-consistent mean-field (SCF) model, which allowed calculation of the polymer concentration profile in the direction normal (z) to the sorbent surface. We will refer to the original SCF model with concentration variation in one dimension as the ID model. An extension of the Scheutjens and Fleer model by Cosgrove et al. [1987a] gives the density distribution of end-grafted chains away from the grafting surface (normal to the interface). This work was shown to give very good quantitative agreement with experimentally measured density profiles [Cosgrove and Ryan, 1990b]. Leermakers et al. [1990] extended this model to a cylindrical lattice (Figure 2.2 (b)) to allow calculation of the segment density distribution in two dimensions (2D model). This allowed Leermakers et al. [1990] to study the energetics of membrane formation and the distribution of chains in inhomogeneous membranes. Our goal is to model a particle interacting with end-grafted polymer chains, which requires both of the above extensions to the original ID self-consistent-field model of Scheutjens and Fleer. In particular, the cylindrical coordinate system of Leermakers et al. [1990] allows the distribution of end-grafted polymer around a particle to be monitored in the direction normal (z) 38 to the grafting surface and in the radial (r) direction (parallel to the grafting surface). In addition, by combining step-weighted walk statistics with a statistical mechanical treatment of the lattice, important thermodynamic information about the system can be calculated such as how much energy it takes to move the particle into the grafted polymer. t r Z v Figure 2.2: Cartesian (a) and cylindrical (b) lattice geometries used in the self- consistent-field model. Areas treated with mean-field averaging are shown in grey as well as the grafting surface in (b). The derivation of our model equations follows that for the original ID model [Fleer et al., 1993; pp. 171-180]. However, due to the addition of an extra dimension there are some important differences, which will be outlined in the following sections. Moreover, all model equations related to the movement of a particle into the grafted chain volume are new. The grand canonical partition function S can be written as a summation of canonical partition functions Q weighted with the appropriate Boltzmann factors: N/i sU,F,r) = £ e ( A ^ , 7 V r (2.11) N 39 The equilibrium distribution of polymer corresponds to the maximum term of Equation (2.11), taken with respect to the number of chains in various conformation distributions {«,c} . This is found from the total differential where n. is the number of molecules of component / adopting a conformation denoted by c. The lattice and polymer conformations in it are further constrained by the requirement that all lattice sites be occupied, which is given by Z Z fc<(^-MO} =0 (2.13) i c where rtc(z,r) is the number of segments of the chain of type / in conformation c that are in the radial of the cylinder in position (z,r). L(r) is the number of lattice sites in radial r. The maximum of a constrained function can be found using the method of Lagrange multipliers. A new function is therefore written / = lnH-ZE«(z^)(SE«(^'-)-^)l <2-14) z r y i c J where afar) represents the Lagrange multiplier for each mean-field radial position in the lattice. Taking the derivative of this expression with respect to one molecule of type k in conformation d gives: 51nS - I Z « M ' / M = 0 (2.15) 40 The grand partition function E is then expanded into its individual terms and evaluated with respect to a reference state denoted by *. The most convenient reference state is that of pure, expansion gives: 0 (2.16) Here Q represents the multiplicity, or the number of ways that a specific set of molecules with fixed internal energy U can be placed on the lattice, and ju is the chemical potential. It is now necessary to develop expressions for Q/Q*, U-lf and JU-JU* based on the lattice model. 2.2.1 Cylindrical Lattice Geometry In the cylindrical lattice (Figure 2.2 (b)), cylinder axis layers are numbered 2=1, ...,Zmax in the direction normal to the grafting surface and rings in the radial direction are numbered r = 1, ...,Rmax- In each ring, within a given layer z, the lattice sites are considered to be indistinguishable and treated by mean-field averaging. The grafting surface at z = 0 is treated as impenetrable and all other surfaces are treated as reflecting boundaries. In all calculations, the lattice is large enough that boundary effects on the chains are negligible. The area of a given radial, specified in lattice units, is given by: unmixed, amorphous components, where S*=J~jH*, Q* = ]~[fi* and U* = ̂ U*. The (2.17) The circumference S(f) is then given by: S(r) = 2nr (2.18) 41 Evaluation of the step-weighted random-walk chain configurations requires knowledge of the step probability for a move from a given lattice site to any adjacent site. We follow the same rules as Leermakers et al. [1990], and allow steps up (higher z), steps down or steps in the same layer. We also allow steps inward (lower r), steps outward, or steps in the same radial shell, and any combination of these two steps simultaneously (i.e. in and down). Physically the lattice step probabilities are a consequence of the assumed geometry and are thus constants. In the normal direction, if a hexagonal lattice is assumed (12 nearest neighbors), the probability to step up or down a layer is X\(z) = X-i(z) =3/12, corresponding to the number of nearest neighbors in the layers adjacent to layer z. The probability to stay in the same layer z is then ô(z) = 6/12. For steps inward and outward the probabilities are made proportional to the lattice site surface area on the face that is crossed so that: A_1(r) = -S(r-l)/L(r) 4 (2.19) \(r) = U(r)lL(r) Thus, the probability to step down and outward simultaneously is given by (2.20) and likewise for all other step directions. 42 2.2.2 Combinatorial Entropy The multiplicity Q. can be calculated from the number of ways of placing volumeless chains in the lattice corrected in a mean-field approximation to account for chain volume. The number of ways o,c of placing a single volumeless chain i with conformation c (which may include a number of degenerate states) is given by < = A ° f t Z A ° f c - W . - V i ) (2-21) where Z is the lattice coordination number, Lct is the number of sites available for placement of the first segment of the chain, and Xc (zs - zs_x, rs - rs_x) is the step probability for insertion of the next segment in the chain. The product is over all segment positions, denoted by a ranking number s up to the full chain length TV,. The number of ways of placing all nf chains can be generalized in the following way to account for all components in all possible conformations: «'=nn^ ( 2 2 2 ) The numerator of Equation (2.22) accounts for the ways of placing all chains of type / in conformation c, and the denominator corrects for the fact that each chain of type ;' is indistinguishable. What remains is to reintroduce chain volume into the "volumeless" multiplicity function Q' defined by Equation (2.22). This is achieved by multiplying O' with the fraction of vacant sites available to each new segment added to the chain. The correction factor for the first segment placed in radial (z,r) is L(r)/L(r) = 1. For the second segment placed in the same 43 radial, the correction factor is (L(r)-\)/L(r). The total correction for filling all the lattice sites in a given radial (z,r) is then ^yT, x i ( r ) • Thus, the total correction to a for occupancy of all possible radials Rmax and all layers Zmax is: nrt«% r / w-n Uf) (2.23) Combination of Equations (2.22) and (2.23), which is our next step, can be greatly simplified by subtracting the multiplicity of the reference state (the result being the total increase in the number of possible molecular arrangements due to mixing). According to Flory [1953], the multiplicity of the pure amorphous polymer chain (see Equation (2.5)) is given by: ._0v>,)!f Z ^ a = (2.24) Combining Equations (2.22) to (2.24) and application of Stirling's approximation then leads to In Nrf (2.25) after careful rearrangement of terms. Here Ac is the product of step probabilities for placing all of the segments in a chain defined by the conformation c. This expression represents the entropic part of the partition function. Taking the derivative of Equation (2.25) with respect to ndk gives the following expression a dnt In \Nkndk j -1 (2.26) which provides the first term of Equation (2.16). 44 2.2.3 Internal Energy In this work, we assume that all energetic interactions are adequately described by the Flory interaction parameter Xy- Thus, the change in energy, with respect to the reference state, for one segment of /' to be placed in radial (z,r) is given by ^ff^Z^jM) (2.27) where the summation on j excludes component /, and the angular brackets represent a layer average. The layer average is needed in SCF theory because all lattice radials adjacent to a given radial (z,r) contain different volume fractions of each component. Thus, the result differs from Flory-Huggins theory, in which the whole solution is assumed to have a uniform distribution of each segment type. The layer average is given by (O }.(z,rj) = S I V ^ ^ O (2-28) z' r' where the summations are over all lattice positions adjacent to radial (z,r), and Xz,_zr._r is the fraction of surface area a site in radial (z,r) has in contact with a site in radial (z',r'). The total energy expression (U-lf)lkT can be obtained by multiplying Equation (2.27) with the number of segments of /' in a lattice radial L(r)^>t(z,r) and summing over all components U-U* 1 kT 2 Z Z Z ( R ) E X * « ( * » R ) * * ( * A R ) ) (2-29) z r where the factor of Vi is needed to account for double counting of pair interactions. An adsorbing interface can be accounted for in the model by including the surface in the summation 45 over all components. Adsorption then occurs if a segment is in a layer next to a surface, such as a segment in layer 1 next to an impenetrable surface at z=0. The surface is treated as a component with 0=1 in layer z-0. Finally, the derivative of Equation (2.29) with respect to ndk can be taken in order to generate the second term of Equation (2.16). The result is: dlu-uy) 1 A Z L J = Z Z Z ^ ^ ^ ^ . ( ^ 0 > (2-30) 2.2.4 Equilibrium Segment Density Distribution In order to complete Equation (2.16), we need an expression for the chemical potential of component k. Evers et al. [1990] have noted that jUk for a multicomponent mixture is most easily derived using the relation A-A* = - A n n Q/Q* + (U - U*) . When this is applied to the bulk solution, we find = 1 + ln** - N k ^ - ^ T L ^ X ^ +NkXz*V (2-31) where the superscript b denotes a bulk volume fraction. This result may also be obtained using the Flory-Huggins expression generalized to ; components, since the more complex combinatorial and energy terms derived above reduce to those of Flory-Huggins theory when applied to a homogeneous bulk solution. Substituting Equations (2.26), (2.30) and (2.31) into Equation (2.16) allows one to write an expression that satisfies the condition of equilibrium 46 In ' ^ ] - Z Z ' / ( ^ % ^ = 0 (2.32) V NX J z r kT where and w,Yz,r), defined by Scheutjens and Fleer [1979,1980], is the mean segment potential, which consists of two separate contributions: uk(z,r) = u'(z,r) + uknt(z,r) (2.34) The first term in Equation (2.34) represents a hard core potential that is independent of segment type and the second term defines an interaction potential that depends on the energetic interactions between segments in the lattice. By careful rearrangement of Equation (2.32), we can determine the number of chains i in a conformation c *" i z r where Gj(z,r) is a Boltzmann factor given by: 47 G,(z,r) = exp( M ' ( r '% r) (2.38) G,(z,r) is also called the free segment weighting factor and represents the preference a free segment of component /' has for being in radial (z,r) with respect to the bulk solution. The free segment weighting factors G, and the step probabilities X allow one to calculate the statistical weight of any conformation. In practice, it is not efficient to generate the weighting for individual conformations since it does not provide an average distribution over all possible conformations. Scheutjens and Fleer [1979] used a matrix formalism developed by Dimarzio and Rubin [1971] to generate simultaneously the statistical weighting of all possible conformations. This method has taken the form of a generalized recurrence relation [Evers, 1990; Leermakers, 1990] that gives the chain- end distribution function as G/(z,r,5 + l|l)=G,(z,r)(G,.(z,r,5|l)) (2.39) where G,(z,r,5+l|l) is the statistical weighting factor (chain end weighting factor) for a chain of s+1 segments to end in radial (z,r) after starting from an end segment with ranking 1. For the case of grafted chains, the chain-end distribution function must be generated from each end of the chain (fixed and free) since inversion symmetry does not apply to end-grafted chains [Cosgrove et al., 1987a]. 2.2.5 Chain Grafting Grafting chains to the end of the lattice at layer 1, is achieved by applying the recurrence formula in Equation (2.39), started from segment 1 with the restriction that: 48 Gi(z,r,\\\) = Gi(z,r) ;z=l,allr' ' (2.40) G,(z,r,l|l) = 0 ; all other z and r (2.41) Grafting positions are specified through different restrictions on Equation (2.40). It is possible, for example, to graft a single chain to the center radial of the cylinder by specifying the only nonzero starting value in Equation (2.40) at (z=l,r=l). The complementary chain-end distribution function is then generated (starting at the other end of the chain) with no restrictions on the location of the free end. This is done using Equation (2.39) from the starting point Gt(z,r,N\N) = G,(z,r) for all z and r. Scheutjens and Fleer showed that the statistical weight of an inner chain segment can be generated using the chain-end weighting factors for a shorter chain of length s and a second chain of length N - s. This approach follows from the "composition law", which is shown schematically in Figure 2.3. s segments N-s+1 segments Figure 2.3: The probability that segment s ends in radial (r,z) is found by combining the probability of two shorter chains ending in that radial. This involves combining two step-weighted walks, one for 1,2,...s segments and one for N,N-l,...s. To calculate the volume fraction ®i(z,r) of a component i in radial (z,r) in our model, the chain-end weighting factors for all chain conformations are computed and combined according to the composition law to give: 49 , ^ v^G,(z,^sl)G,(z,r,.s^ % & r) = C, 2̂ ^rr—, — (2.42) The summation s is over all segment positions in the chain and a normalization constant Q is needed to transform the weighting factors to volume fraction. The normalization constant for free chains is then given by Q,= <D* / N,,. For components, which are free to move in or out of the system, <D* must be specified. If however, one component is restricted to remain within a certain area of the lattice, such as end- grafted chains, the true bulk volume fraction O 6 will be zero for that component. It then becomes necessary to have a different way of specifying the amount of that component in the system. Scheutjens and Fleer [1985] therefore introduced the idea of restricted equilibrium. Under restricted equilibrium conditions, one can replace <D6 with a pseudo-equilibrium bulk volume fraction calculated using the chain end weighting factors. Summing the equilibrium expression (Equation (2.37)) over all possible conformations «,c gives «,. Substitution of Q = Of I Nt then allows one to write the normalization constant as: C , = j:^L(")G!(z,r,N\l) <"3> This expression can be used with Equation (2.42) for molecules in restricted equilibrium, where rij is the number of molecules of/, which must be specified for the calculation. The denominator of the expression gives the total statistical weight to find a set of chains containing all possible conformations in the cylindrical lattice. 50 The combination of Equations (2.38), (2.34), (2.36) and (2.42) constitutes a set of self- consistent equations that must be solved using numerical techniques. This is because the volume fraction 0;(z,r) of each component is calculated from the free segment weighting factors Gt(z,r) based on the mean segment potentials Ut(z,r). The mean segment potentials in turn depend on the volume fraction distribution of the components through the functional dependence of the interaction potential Ujint(z,r). 2.2.6 Interaction Energy The partition function Q is related to the change in Helmholtz energy of the system to by ln(e/e-) = % ^ = l n ( n / n - ) - M (2.44) where the second equality separates the Helmholtz energy into its energetic and entropic component: S-S* =-k\n(Cl/Q*) (2.45) Equation (2.45) allows one to define the mixing entropy in terms of the free segment weighting factors and volume fractions by substitution of the equilibrium expression (Equation (2.37)) into the combinatorial expression (Equation (2.25)). Taking the summation over all conformations gives: {s-s^/k^ZT^^M 1 2 r ln<D* —-i-+lnG,(z,r) (2.46) By adding the energy term to this expression the complete expression for the change in Helmholtz energy is obtained: 51 z r i ln<D' UlnG,(z,r)+i2^(oy(z,r)) (2.47) The restricted equilibrium criteria then allow one to make a useful substitution in Equation (2.47) for the bulk volume fraction of grafted components: Finally the energy change due to the movement of a particle into the grafted polymer chains is calculated by Am\z) = \A-A*\z)-{A-A*\oo) (2.49) where the coordinate z represents the layer position of the leading edge of the particle. The particle is moved in the z-direction toward the surface one layer at a time and the brush distribution and Helmholtz energy are calculated at each position. 2.2.7 Numerical Solution The set of self-consistent equations may be solved numerically by writing a set of non- linear functions [Leermakers, 1988] of the form: , 2 > , ' ( ^ ) i i The set includes an equation for each component in each mean-field radial of the lattice. The first two terms of Equation (2.50) require that the volume fractions of all components add up to unity in each lattice site, and the second two terms are only satisfied when the hard core potential 52 becomes independent of segment type. An example of the computer code written for the calculations is provided in Appendix A. The problem could be solved using a global minimization method, where the root of the sum of squares F = f,2(z,r) z r i 1/2 is minimized, but a non-linear equation method has been chosen, since the functions are fairly well behaved, and the need for second-order derivatives is avoided. No conditions have been observed where the Jacobian / ft ft ft, J = dux ' du2 ' du3 ' ft2 ft2 ft2 (2.51) becomes singular (i.e. local minima) making it possible to carry out iterations until the sum of squares is below a specified tolerance of (typically) le-14. This generally results in volume fractions that are correct to 5 or 6 significant figures. Because of the large number of equations generated when even a small lattice is used, it is a challenge to devise a rapid solution scheme (i.e. 25 layers by 10 radials with 3 components gives 750 variables). It is common to solve a set of non-linear equations using Newton's method, or some form of modified Newton's method [Scales, 1985]. Whatever method is used, it is necessary to have some gradient information about the functions being solved. This is a very costly step if the gradient is evaluated numerically each iteration, since finite differences must be carried out for each function with respect to each iteration variable. In the class of 53 methods known as Quasi-Newton methods, the information obtained from a previous iteration is used to update the Jacobian (or Hessian) matrix so that the solution is more rapid. In a typical Newton iteration, it is necessary to solve a large set of linear equations in order to calculate the search direction. This can sometimes be avoided when the Jacobian is tri- diagonal or has some specific symmetry. Decompositions might also be used to speed up the solution of the linear equations. With a large set of equations, the solution by Gaussian reduction can be the single most time consuming step. It is therefore often desirable to approximate the Jacobian for the first iteration, invert the matrix and obtain the search direction by direct multiplication. If an updated inverse is then calculated, the method will proceed rapidly toward the solution, even for the case of several hundred variables. An algorithm having these time saving features is known as Broyden's rank-one method, and uses an update formula developed to improve the inverse Jacobian [Broyden, 1965]. An outline of the algorithm is as follows [Scales, 1985; p. 128]: 1. Input u, TOL 2. Approximate J, and calculate J ' 1 3. Calculate An = - J 1 *f 4. Compute a {so that F(u + a* Au) < F(u) } 5. Setu = u + a*Au 6. Compute Q J 7. Set/'1 = J A + Q J 8. Repeat steps 3 to 7 until F < TOL In this algorithm, J and J ' 1 are approximations to the Jacobian and its inverse respectively. The calculation of a is done using a line search algorithm, in order to minimize the sum of squares 54 (approximately) in the search direction, and the matrix Q J is calculated using Broyden's update formula: ( / - ' A / - A K ) A K ^ * AuTJ'Af K } Broyden's update formula comes from an extension of the secant method to multiple dimensions. In a one dimensional secant method, the finite difference approximation of the first derivative is replaced by the difference in the last two function evaluations over the difference in the last two values of the independent variable: ( 2 5 3 ) This saves time if the extra function evaluations needed for finite differences are costly. The method can be extended to multiple dimensions, where the analogous model becomes: /(« t) = / K - . ) + ̂ i K - « H ) (2-54) In multiple dimensions Equation (2.54) does not provide enough information to specify / precisely, because as long as AM is not equal to zero, a large number of matrices / can satisfy Equation (2.54) [Dennis, 1983]. Broyden's method, however, generates J (the Jacobian approximation) by minimizing the change in the model while still satisfying Equation (2.54). Broyden's [1965] Jacobian update formula is given by k+] k AuTAu which, when rearranged using Householder's formula [Dennis, 1983], gives the inverse update formula shown as Equation (2.52). 55 The search direction that is generated using Newton's method may not always reduce the sum of squares if the full step is taken. When attempting to solve a set of nonlinear equations by minimizing their sum of squares, it can be shown that the Newton direction should always produce a descent direction [Dennis, 1983], as long as a local minimum is not encountered. For this reason it is often desirable to take a fractional step in the search direction. The subroutine LINSRCH carries out a univariate minimization, in which the sum of squares is minimized with respect to the fraction, a, using the current search direction. Brent [1973] first published the algorithm. The minimization interval of the function is specified to be between XI and XF assuming a unimodal function. Brent's [1973] subroutine uses a hybrid method, which combines a quadratic interpolation with "Golden Section" search. If the predicted position of the next point found by quadratic interpolation does not meet several criteria, such as falling between the last two lowest points, being non-zero, and not being too small or large a step, then a "Golden Section" step is taken. This is a fractional step toward the lower of the points. By combining the two methods, the algorithm is guaranteed to find the minimum of a function as long as one exists on the search interval. In practice, the fraction a is chosen to lie between 0.001 and 0.7 and the search is allowed to proceed for a maximum number of interpolations, usually 5 to 7 before accepting the change in u of a*Au. There is also a termination criterion based on the width of the search interval (which narrows as it proceeds) and the distance of the next interpolation from the interval edges. Our algorithm is somewhat simplistic, which warrants some discussion given the complex problem to which it has been applied. The use of a line search helps to make the 56 solution method more global. However, convergence is not assured when a poor initial guess of the starting values is made. An increase in the sum of squares can be observed in cases where a poor search direction is generated by the Jacobian update, since the line search terminates after a fixed number of iterations. In general though, the solution will wander uphill slightly under such conditions and then make large steps downhill from the new position. 2.3 Evaluation of Model Parameters 2.3.1 Scaling the Lattice Model It can be difficult to describe the results of lattice-model calculations in terms of real dimensions. This is due to the fact that a real molecule, with many degrees of conformational freedom and non-uniform shape, is being described by a set of contiguous lattice sites, which have a fixed volume and simplistic spatial arrangement. When a solution of different types of molecules is being described by a single lattice geometry, the lattice scaling is made even more difficult by the desire to find a lattice size that simultaneously describes well the different types of molecules. In this study, the physical properties of the system that require scaling include the polymer chain length, the surface density of polymer and the size of the particle (protein) interacting with the polymer. The lattice analogs of these properties are summarized in Table 2.1. There are different methods of scaling, depending on the lattice model chosen. In the original work by Flory dealing with a bulk solution of polymer, the dimensions of a lattice unit were never required. Polymer chain length in the lattice was scaled against the size of the solvent molecules. The number of segments in a polymer chain was then calculated as the ratio of the polymer to solvent molar volume [Flory, 1953]. In more complex lattice models, such as 57 ours, where concentration gradients exist or chains are tethered to a specific location, more complicated scaling techniques are necessary. Table 2.1: Experimental properties that require direct comparison with lattice values. Property Lattice Value Experimental Value Chain Length Surface Density Particle Size N (segments) 0 (monolayers) R (lattice units) r (monomers) r(mg/m2) **(A) 2.3.1.1 Lattice Unit and Chain Length Meaningful comparison of a lattice chain to a real chain is possible by making use of two chain properties, the radius of gyration and the contour length. Equating these two properties between lattice chains and real chains allows one to calculate the appropriate dimension for one lattice unit and the number of lattice segments N corresponding to a chain of a given molecular weight. In experimental systems, the radius of gyration is given by Rg=aMU2=a(rMmf2 (2.56) where a is an experimentally determined parameter and M is the polymer molecular weight, also given by the number of monomers r multiplied by the monomer molecular weight Mm [Fleer etal., 1993]. For chains described using a random-walk model, the radius of gyration is given by Rg=pNBl2 (2.57) 58 where NB is the number of bonds per chain (NB = r -1), / is the bond length and p is the persistence (a measure of the chain stiffness). If we treat the smallest step in the random-walk as a monomer, then equating Equations (2.56) and (2.57) give a simple expression for p based on experiment a2rMm a2Mm „. P = T~ « —T2- (2.58) F NJ2 I2 V ' since at large NB, r INB «1 . The persistence pL can also be calculated directly by the for lattice model 1 P l = 6 1-z-1 (2.59) where z is the lattice coordination number [Fleer et al., 1993]. Here pL is a consequence of the lattice geometry only. Equating the random-walk radius of gyration (Equation (2.57)) for a real chain with that for a lattice chain gives pNBl2 = pLNLBl2L (2.60) where TV^ is the number of lattice bonds ( = N -1) and lL is the segment length, or length of one lattice unit. Similarly, equating the fully extended chain length (contour length) of both chains gives: NBl = NLBlL (2.61) 59 Combining Equations (2.60) and (2.61) and substituting for the persistence values (Equations (2.58) and (2.59)) of each chain then allows us to calculate the dimensions of the lattice unit based on experimental parameters (2.62) PL Finally, substituting this result back into Equation (2.61) allows us to calculate the lattice chain length N. For example, polyethylene has an Mm= 28, a= 0.434 nm [Fleer et al., 1993], and /= 0.25 nm; in a cubic lattice, pL= 0.233. This gives lL- 0.9 nm and a chain consisting of 100 monomers should have a number of lattice segments N equal to about 28. In cases where the monomers are large, or the chains are very stiff, this method may break down. The resulting lattice unit sizes can become too large to allow simultaneous scaling of the solvent size and polymer size. A second method of scaling the chain length is based on rotational isomeric state theory [Flory, 1969]. Cosgrove and Ryan [1990b] have scaled their results based on work by Cohen Stuart [1980], which relates NB to Nw by: NB/NLB= log z (2.63) Here z is the lattice coordination number and zp is a characteristic constant related to the chain persistence. It is also possible to arbitrarily specify the dimension of a lattice unit. Baskir et al. [1989a; 1989b; 1987] chose a lattice unit of 4A to model polyethylene oxide and Dextran, and 60 scaled the number of chain segments equal to the degree of polymerization. Using this simple approach, they obtained good quantitative comparison with aqueous two-phase partitioning data for various globular proteins. 2.3.1.2 Polymer Surface Concentration Once a lattice unit is chosen and a lattice chain length defined, the polymer surface density can be calculated. However, due to the limitations imposed by treating real chains as segments in an ideal lattice, comparing lattice graft densities to experiment is not trivial. In the simplest case, a lattice chain can be considered to occupy a real volume equivalent to that of the lattice sites it occupies. In this case, the surface density Y (mg/m2) is related to the lattice analog 6 by where m (mg/chain) is the mass of one polymer chain. In comparison to experiment, however, Equation (2.64) is not always satisfied. For example, if a lattice unit of 1.0 nm is used (100 A2 per surface site) and it is experimentally determined that the chains are grafted with 100 A2 per chain, the lattice model requires that the grafting density <r be unity, indicating no solvent is present in the lattice model. Clearly this is unrealistic, so at high surface densities it may be necessary to arbitrarily reduce the value of lL in order to reduce <J . If lL is made too small, however, the latice chains can have an unrealistic density. This occurs because, although the contour length is matched, the cross-sectional area of the lattice chains is not scaled properly. The problem of matching lattice chains exactly to the volume of real chains can be avoided by scaling the total number of grafted segments based on constant polymer density. 61 Cosgrove and Ryan [1990b] scaled their lattice model in order to compare calculated polymer density distributions of end-grafted polyethylene oxide (PEO) with those measured using small- angle neutron scattering (SANS). They used the true polymer density ppolymer and the lattice unit to calculate a monolayer coverage mon = ^LP polymer (2-65) in units of (mg/m2). It then follows that the experimental coverage is equal to the monolayer capacity multiplied by the number of equivalent monolayers #of segments in the lattice . r = KorP (2-66) Cosgrove and Ryan's method gave very good quantitative agreement between experiment and model calculations. They state, however, that the method may break down at high surface overage, or for very high molecular weight chains. 2.3.1.3 Particle Size In order to model a particle interacting with the end-grafted polymer, a third component has been added to the lattice with a specific geometry. In the cylindrical lattice the most convenient particle geometry is that of a small cylinder (usually with length set equal to radius) in which the volume fraction of particle is specified to be very high ( ® p > 0.99). The radius of the cylinder then allows scaling to the relative size of real molecules, such as a protein with a specific molecular weight. Thus, we assume that a particle with a given radius has about the same interaction with the polymer as a protein of the same radius would. This seems reasonable because, in terms of particle dimensions, our calculations show that the radius has the dominant effect on brush-particle interaction energies. 62 Tyn and Gusek [1990] cite several different empirical correlations for estimating protein diffusion coefficients using protein molecular weight or radius. By coupling two of these methods it is possible to derive an equation, which directly relates molecular weight to radius. For example, based on the Stokes-Einstein equation, Poison developed a correlation relating the diffusion coefficient of a protein to its molecular weight D = A/MVi (2.67) where A is a constant and M is the protein molecular weight. Tyn and Gusek proposed the following correlation between D and the radius of gyration Rg of the protein D = B I Rg (2.68) where B is a constant. Equating these two expressions gives M = {CRg)3 (2.69) where C is a constant that can be fit over a molecular weight range of interest. For proteins up to about M = 150 kDa, a non-linear curve fit was carried out using Equation (2.69). Figure 2.4 shows that a very good fit was obtained, giving a value of C = 1.50 + 0.02 (kDa1/3/A), allowing direct scaling of model particles to real molecular weights. 63 10 15 20 25 30 35 Radius (A) Figure 2.4: Non-linear curve fit of data taken from Tyn and Gusek [1990] for use scaling particle size to protein molecular weight. 2.3.2 Measuring Interaction Parameters The Flory interaction parameter xtj is based on the statistical mechanical quantity known as the interchange energy Acoy. The interchange energy is a physico-chemical property, which can and has been measured experimentally. Several techniques can be used to measure chemical potentials of solvent and/or polymer including differential-vapour-pressure, osmometry and low- angle laser-light-scattering (LALLS) [Rathbone et al., 1990; Haynes et al, 1989]. Based on predictions of the chemical potential using our model, interaction parameters may then be calculated from the experimental results. 64 2.3.2.1 Differential-Vapour-Pressure and Osmometry To fix ideas, consider a binary mixture of solvent (1) and polymer (2). For such a mixture, the derivative of the Gibbs' energy of mixing from Flory-Huggins theory taken with respect to the number of solvent molecules [dGm I dnx \ p n gives the solvent chemical potential: Mi ~M° =RT ln(l-<D2) + <D + Xl2®2 (2.70) The chemical potential of a component is related to its activity by fix -ix° = RTlna] (2.71) where ax is the activity of the solvent. In the differential-vapour-pressure experiment, a, is related to the difference in the vapour pressure between pure solvent and the polymer-containing solvent A P by (Pi°-&P) , , 1 where P° is the vapour pressure of pure solvent. If P° is low, the vapour phase can be treated as ideal. Haynes et al. [1989] give a detailed description of the experimental technique and apparatus used for differential-vapour-pressure measurement. Substitution of Equations (2.70) and (2.72) into Equation (2.71) with subsequent rearrangement allows direct calculation of X12 for a given concentration of polymer. It is useful to note that the polymer volume fraction $ 2 is equal to the weight concentration C2 (g/mL) multiplied by the partial specific volume of the polymer in solution v2 (mL/g). 65 Osmotic pressure measurements can be used in the same manner as vapour pressure measurements to calculate interaction parameters. The osmotic pressure of the solvent (2.73) is related directly to the solvent chemical potential and the molar volume of the solvent Vi. 2.3.2.2 Low-Angle Laser-Light-Scattering ( L A L L S ) A more direct measure of nonideal intermolecular interaction in solution can be gained from the light-scattering properties of its molecules. Consider once again a dilute binary polymer solution of solvent (1) and polymer (2). The osmotic virial expansion of McMillan and Mayer [1945] expresses the solvent chemical potential JUJ in such a solution exactly in terms of an expansion in polymer concentration C2 (g/mL) (2.74) where Vl is the molar volume of the solvent, MW2 is the weight-averaged molecular weight of the polymer, and A22 (mL-mol/g2) and A222 (mL2-mol/g3) are the osmotic second and third virial coefficients, respectively. Osmotic virial coefficients are dependent only on temperature and the nature of the solute and solvent at normal pressures. Statistical mechanics suggests that A22 accounts for two-body interactions in solution and A222 accounts for three-body interactions. At low enough concentrations it is possible to truncate the virial expansion after the second-order terms since higher-order interactions will be approximately negligible. ^-^=-RTVx M •-\-A22c2 ^ A222c2 +. W2 66 For dilute polymer solutions, Equation (2.74) can be truncated after the second virial coefficient term without loss in accuracy. The solvent chemical potential jui can then be determined as a function of polymer concentration if MW2 and A22 are known. The values of both of these constants can be determined by low-angle laser-light scattering (LALLS), which measures the light-scattering properties of a solute (polymer) molecule whose size is small compared to the wavelength of the incident light (i.e., Rayleight scattering). A detailed description of the technique and associated theory, which is not of direct relevance to this thesis, is provided by Rathbone et al. [1990]. Equating the right-hand side of Equation (2.74), truncated after the second virial coefficient term, with that of Equation (2.70), allows regression of Xn from LALLS data. One advantage of the LALLS technique over the two colligative methods (vapour-pressure and membrane osmometry) is that it is easily extended to multicomponent measurements, allowing determination of Xtj parameters characterizing the interaction of solute /' with solute j (i.e., polymer /' with protein j). The theory used to calculate x<j parameters from LALLS data for a ternary aqueous polymer solution is that derived in Haynes et al. [1993]. The principle equations relating the osmotic second virial coefficients to the Flory interaction parameters are (2.75) (2.76) where v, is the partial specific volume of component / in solution. 67 Regrettably, Xj data are extremely rare for aqueous polymer systems. LALLS data (as well as vapour-pressure and membrane osmometry data to a lesser extent) are available for a few aqueous solutions containing polymer and protein [Haynes et al., 1993; Rathbone et al., 1990; King et al., 1988]. Equations (2.70) to (2.76) are therefore valuable because they allow determination of the realistic range of Xij values we would expect to observe in our systems of interest; namely, the interaction of a protein with a grafted-polymer brush, both of which are solvated by water. Flory interaction parameters are shown in Table 2.2 to Table 2.6 for a number of polymers and proteins regressed from osmotic second virial coefficients measured by LALLS [Haynes et al., 1993]. All polymer-solvent and protein-solvent interaction parameters (Table 2.1 to Table 2.5) are significantly positive indicating an unfavourable enthalpic interaction between the solute and water. In a very good solvent, X12 is negative or very close to zero (athermal). In cases where xn is positive, solvation can only result because of a large increase in entropy on mixing. Most water-soluble polymers have X12 values close to the theoretical point at which phase separation will occur for polymer of infinite molecular weight {xn ~ 0.5). Thus, decreased solvent quality (increased X'2) in these cases can lead to phase separation of polymer or precipitation of proteins. Interaction parameters X23 shown in Table 2.6 indicate that there is a very weak enthalpic interaction between the polymers (PEG or dextran) and proteins. Bovine serum albumin and lysozyme appear to have a very weak attraction to the polymers while a- Chymotrypsin has a very weak repulsion. 68 Table 2.2: Polymer-solvent interaction parameters for PEG and Dextran in water at 25 °C calculated from light-scattering data of Haynes et al. [1993]. Polymer M w (g/mol) A N X 1 0 4 (mL-mol/g2) X12 PEG 3350 3860 36.3 0.45 PEG 8000 11700 30.3 0.45 Dextran T-70 10100 4.0 0.49 Dextran T-500 43400 1.3 0.50 Table 2.3: Protein-solvent interaction parameters calculated from light-scattering data for bovine serum albumin (BSA) in water at 25 °C. Buffer Salt (or Salt) Salt Cone, (mol/kg) Ionic Strength (mol/kg) PH A N X 1 0 4 (mL-mol/g2) 30 2 Potassium Phosphate — 0.075 7.0 8.3 0.47 — 0.14 7.0 6.1 0.48 — 0.28 7,0 3.7 0.49 — 0.4 7.0 2.5 0.49 — 0.5 7.0 2.4 0.49 KCI 0.05 — 7 8.7 0.47 NaCl 0.05 — 7 8.9 0.47 Table 2.4: Protein-solvent interaction parameters calculated from light-scattering data for lysozyme in water at 25 °C. Buffer Salt Salt Cone. Ionic Strength A H X 10 4 X12 (or Salt) (mol/kg) (mol/kg) PH (mL-mol/g2) Potassium Phosphate — 0.075 7.0 4.1 0.49 — 0.14 7.0 -2.5 0.51 — 0.28 7.0 -7.2 0.53 — 0.4 7.0 -10.0 0.54 — 0.5 7.0 -10.0 0.54 KCI 0.05 — 7 5.8 0.48 NaCl 0.05 — 7 6.1 0.48 69 Table 2.5: Protein-solvent interaction parameters calculated from light-scattering data for a-chymotrypsin in water at 25 °C. Buffer Salt Salt Cone. Ionic Strength A N X 1 0 4 X i 2 (or Salt) (mol/kg) (mol/kg) pH (mL-mol/g 2) Potassium Phosphate — 0.075 7.0 -26.0 0.59 — 0.14 7.0 -13.0 0.55 — 0.28 7.0 -4.9 0.52 — 0.4 7.0 -3.5 0.51 — 0.5 7.0 -2.8 0.51 KCI 0.05 — 7 -13.0 0.55 NaCl 0.05 — .7 -13.0 0.55 Table 2.6: Protein-polymer interaction parameters X 2 3 calculated from light-scattering data in potassium phosphate buffer at 25 °C. Protein Polymer P E G 3350 P E G 8000 Dextran T-70 Dextran T-500 BSA pH7 : 50mM buffer pH7 : 100mM buffer -0.07 -0.06 -0.08 -0.06 -0.04 -0.03 -0.04 -0.02 Lysozyme pH7 : 50mM buffer pH7 : 100mM buffer -0.06 -0.04 -0.06 -0.04 -0.02 0.00 -0.02 0.00 a-Chymotrypsin pH7 : 50mM buffer pH7 : 100mM buffer 0.06 0.03 0.06 0.03 0.10 0.07 0.10 0.07 70 3 Isolated Terminally-Attached Chains 3.1 Motivation Before proceeding to our analysis of the interaction of a particle (protein) with a brush, it is useful to gain a better understanding of how grafted-chain conformation is altered by the movement of an impenetrable particle or surface into its volume space. The only meaningful study is that of Subramanian et al. [1995], who modeled the compression of a single grafted polymer chain using scaling theory (see Section 3.2). In this chapter, we apply our SCF model to the calculation of grafted-chain configurations and energetics during compression by a cylinder of radius R. A schematic of the system under investigation showing all characteristic dimensions is shown in Figure 3.1. Polymer segment density distributions and compression energies have been calculated in the cylindrical lattice for mushrooms of various chain lengths compressed by discs of varying radius. Jdisc Figure 3.1 ms Schematic diagram of a single end-grafted chain in (a) uncompressed (b) compressed, and (c) escaped states. The chain grafting point indicated by • is at lattice position (z=0,r=0). Characteristic dimensions include: (i) disc position, Zdisc, (ii) disc radius, RdiSC (hi) r.m.s. radius of segments from cylinder axis, Rrms and, (iv) Flory radius, RF. 71 3.2 Scaling Theory of Mushroom Compression Subramanian et al. [1995] model the compression of a single polymer chain by a finite disc of radius Rdisc using scaling theory. They have studied the case where the disc is centered over the grafting point of the chain, and as a result the chain cannot fully escape under compression. Under such conditions, they predict that the chain will undergo pure compression until a point is reached at which a first-order transition occurs, where the free end of the chain suddenly escapes from under the disc. The scaling development is briefly described here. Subramanian et al. [1995] consider the grafted chain to be a mushroom-type coil of radius Rg = IN3'5 which can be modeled as a string of blobs each of diameter D under compression, where D is also the distance of the disc from the grafting surface. The number of blobs B (each with TVB segments) in the chain are calculated from the scaling result that the blob radius is given by D = IN3^5 and thus B = NINB = N(lID)5'3. The formation of each blob is associated with energy kT and the total energy associated with chain compression is thus Ablob=kTN(l/D)5/3. Partial chain escape is predicted by the model at conditions where the chain can lower its energy by stretching a tether to the edge of the disc and forming a larger escaped blob outside the disc radius. If the tether contains Bt blobs, the tether energy based on unperturbed blob formation is kTBt. However, the tether may stretch to allow partial chain escape from the disc radius Rdisc, particularly when Rdisc is larger than the two dimensional radius of the string of blobs in the tether, given by R1DMher = DB3IA. The stretching energy required for this process is given by scaling arguments to be kT(Rdjsc I R t Y • The total free energy of the mushroom may then be written as ATotal = fcT(Bt + (Rdisc/D)4B;3). 72 The minimum total free energy with respect to the number of blobs in the tether is found when B, -3xl4(Rdisc ID). Subramanian et al. [1995] equate this minimum total energy to the energy of the fully confined chain (i.e. no escape or stretching) to determine the critical point for chain escape, given by Dc =(31/2 142/3)N3/2l5'2 IR3,'2. The two dimensional radius of the chain just before escape is therefore predicted to be R2D - DCB3'4 - Rdisc(Rg IRdisc)sn . From this expression, Subramanian et al. predicted that if the disc radius is greater than the free coil radius of gyration, there will be a jump in the chain radius upon escape that is analogous to a first-order transition. They also predict that there is an energy barrier to be overcome for chain escape to occur, which is equal to the energy required for the chain to stretch to the edge of the disc. This stretching energy is given by kT(Rdisc IR2D(DC))4 and can be significantly larger than kT. Using similar scaling arguments Guffond et al. [1997] studied the compression of end- grafted mushrooms in a poor solvent by various flat discs and step-discs finding similar results for an escape transition. Williams and MacKintosh [1995] also studied the compression of mushrooms under curved surfaces and reported regimes of surface curvature resulting in an escape transition and an unconfined chain. They also investigated compression under irregular surfaces finding multiple-escape transitions. 3.3 Analysis of Mushroom Compression using Numerical SCF Theory At low grafting density, the conformation of an end-grafted chain is based only on the impenetrable nature of the grafting surface and interactions of the chain with itself and solvent, assuming no other species are present. As mentioned previously (Chapter 1), de Gennes [1980] has shown that a single end-grafted chain in a good solvent forms a random coil, roughly in the 73 form of a half sphere, with a Flory radius of approximately RF « / J V 3 / 5 . Thus, a chain of 100 segments will have a Flory radius of about 16 lattice units. The root-mean-squared distance of segments from the cylinder axis R~... - 1/2 (3 .1) N has also been calculated in this work in order to characterize the amount of spreading and escape a chain undergoes when compressed. The fractional segment distribution in both the radial (r) and normal (z) directions is shown in Figure 3.2 for an end-grafted chain of 100 segments in a good solvent ( #2=0) under no compression. The Flory radius and Rrms of the coil as calculated in Equation (3 .1) are also shown in Figure 3.2 for comparison. The number of segments in each radial n(r) is calculated as n(r) = Y,L(r)0(z,r) (3.2) z which is simply the volume fraction 3>,(z,r) multiplied by the number of lattice sites L in the radial r summed over all layers. The fractional distribution was then computed by normalizing Equation (3 .2) by the number of segments in the chain. In agreement with de Gennes' predictions, we find a maximum in the segment fraction a few layers from the surface in the normal direction. There is also a maximum in the segment fraction a few radials from the center. This is due to the fact that the inner radials have fewer lattice sites and can thus accommodate fewer segments despite having higher segment densities. The two curves approach zero at roughly the same distance, suggesting the chain forms roughly a half-sphere (or mushroom), spreading more in the radial than normal direction. 74 Layer or Radial Figure 3.2 Fractional segment distribution of a single end-grafted chain in the radial (r) (dash-dot line) and normal (z) (solid line) directions with 7V=100 segments. Dashed line shows the root mean squared (rms) distance of segments from the center axis, RRMS, and dotted line shows the Flory radius, RF. The segment density profile of a single end-grafted chain is shown in Figure 3.3 (a). The density is highest at the point of grafting and drops off rapidly with distance from the surface. The graft density at (1,1) is about 0.4, which corresponds to slightly more than 1 segment in the % lattice sites of the first radial in layer 1. In Figure 3.3 (b), a disc of radius 20 with a thickness of 4 layers has been placed in layer 2 and the chain distribution around the disc is shown. It can be seen that the chain density becomes much higher in the outer radials (r = 15-20) of layer 1. A significant portion of the chain is escaped from under the disc and forms a smaller coil outside the edges of the disc. 75 (a) (b) Figure 3.3 Segment density distribution of a single end-grafted chain with N=100. The chain is shown in (a) with no disc and in (b) with a disc of radius, RdiSC=20, sitting in layer 2 compressing the mushroom. The extent to which segments spread out in the radial direction under a compressing disc is proportional to Rrms (Equation (3.1)). The r.m.s. distance of segments from the center of the cylinder is shown in Figure 3.4 for a chain of 100 segments compressed by two discs of different diameters. In both curves, Rrms is the same under compression and decompression. We see no hysteresis in chain density profiles. The chain distribution is generated in our model using a step-weighted walk starting from position (1,1). If the disc is placed in layer 2 for example, we find only one self-consistent solution for the chain density distribution around that disc. An important feature of Figure 3.4 is that as the disc is made larger, more compression is required in order for chain escape. Consequently, the escape region is sharper due to the fact that more work is required to squeeze the chain to the edge of a larger disc. The escape of the chain under the disc of radius 20 is gradual, as seen by the solid line in Figure 3.4. 76 R. rms R d iso=20 R d isc=30 5 10 15 20 Disc Position (z) Figure 3.4 Root mean squared distance of segments from the center of the cylinder, Rims, is shown as the disc compresses a mushroom with N=T00. Rdjsc=20 (solid), Rdis<=30 (dashed) As discussed in section 3.2, Subramanian et al. [1995; 1996] report that for a disc larger than RF, a first-order escape transition is found. They argue that this transition is associated with an energy required to stretch the compressed polymer chain to the edge of the disc. However, for such a system (a disc of radius 20 compressing a polymer of Flory radius ca. 16), we observe a gradual, reversible transition, which does not conform to a first-order escape transition. The scaling model of Subramanian et al. equates the calculated energy of two limiting cases in order to predict the compression at which a partially-escaped chain is in equilibrium with a compressed chain whose 2D Flory radius is close but not equal to the disc radius. The strongly compressed chain is modeled crudely as a string of blobs, each with energy kT, and the escaped chain is 77 modeled as a strongly stretched string of blobs (tether) with a fraction of the chain escaped. The fact that no such first-order transition is observed in our calculations suggests that information about the transition between these two states is lost in the scaling model, which is not altogether surprising given the assumptions that are made. z disc decreasing 0.10- 10 20 Radial (r) Figure 3.5 Fractional segment distribution in the radial direction as a mushroom with N=75 is compressed by a disc with Rdisc=20. The curves represent the disc in positions of z=16, 8, 5, 4, 3, and 2. The segment fraction in each radial is shown in Figure 3.5 for a chain of A/=75 segments compressed by a disc of radius RdiSC-20. This set of curves shows the results when the disc is placed in layer 16, 8, 5, 4, 3, and 2. Initial chain escape can be seen when the disc is placed in layer 5, and the fraction of segments escaped increases significantly as the disc is moved in further. The Rrms curve for compression shown in Figure 3.6 looks very similar to the dashed line in Figure 3.4, but with a slightly less abrupt escape. It is interesting to note that with the disc 78 in layer 2, the number of covered segments in each radial is about 1.5. Within the mean field approximation, we could interpret this result as formation of a strongly stretched tether under the disc with a large escaped portion. As the chain becomes longer, the tether becomes more strongly stretched and at infinite chain length, would become fully stretched to its limit of one segment per radial under a finite sized disc. . 1 1 1 1 r 5 10 15 Disc Position (z) Figure 3.6 Root mean squared distance of segments from the center of the cylinder, Rims, is shown as the disc of Rdisc=20 compresses a mushroom with N=75. The energetic cost of compressing a mushroom has been calculated based on the loss of conformational entropy. Figure 3.7 shows the compression energy for four different chain lengths. The curves represent compression under a disc of radius i? ĉ=20, and the chain lengths are N=20, 50, 100, and 200 segments. The curves are all monotonically increasing with compression and show no behavior indicating a first-order transition. Comparing the calculated energy curves for long chains with that for the N=20 chain, which cannot escape, shows that 79 chain escape is not accompanied by anomalies in A'"'. The only clear difference between the systems is a shift of the curve to the right (higher disc positions) with increasing chain length, due to the increased excluded volume repulsion of the compressed chain. For the longer chains, JV=50, 7V=100 and JV=200, it can be seen that portions of the energy curves coincide at high compression. The energy curves for N=50 and JV=100 overlap on the 7V=200 curve at a point where the shorter chains have just started to escape. Further compression therefore involves only the work of compressing the polymer tether, which is independent of chain length. Disc Position (z) Figure 3.7 Energy required to compress a mushroom with a disc of radius, RdiSC=20, for several chain lengths: N = 2 0 (dashed), N = 5 0 (dotted), N = 1 0 0 (dot- dash), N = 2 0 0 (solid) In order to investigate the chain escape further, chain length was varied under a fixed disc of radius RdiSC~20. The disc extended from layer 3 to z=Zmax. Figure 3.8 shows the fraction of segments covered by the disc (dotted) and the fraction of segments escaped (solid) from under it. 80 Again, a smooth curve is seen with no discontinuities that would indicate a first-order transition. The onset of escape occurs at just over 50 segments. As a final test for stability of the compressed chain, the chemical potential of the grafted chain was calculated for various chain lengths under a fixed disc and is shown in Figure 3.9 for a chain of JV = 80. Conditions were the same as in Figure 3.8. In the binary mixture under consideration (l=solvent, 2=grafted-chain), phase instability will occur if (5n2/8n2)T,ni,v < 0 [Prausnitz et al, 1986]. However, (8p2/on2)T,ni,v remains positive as shown in Figure 3.10 for all conditions (i.e. chain length and escaped fractions) investigated in Figure 3.8, providing more thermodynamically rigorous proof that although chain escape occurs, it does not involve a first-order phase transition. 50 100 150 Chain Length 200 Figure 3.8 Fraction of segments covered (dotted) and escaped (solid) from under a stationary disc of radius 20, extending from z=3 to the end of the cylinder (Zmax), as chain length is varied. 81 -50.0 -50.5 - -51.0- -51.5- =£" -52.0- -52.5- -53.0 - -53.5 0.4 0.6 0.8 1.0 1.2 1.4 i 1.6 Figure 3.9 Chemical potential of a single end-grafted chain under a stationary disc of radius 20, that extends from z=3 to the end of the cylinder (Zmax), as a function of the number of grafted chains. Figure 3.10 The change in chemical potential with respect to the number of grafted chains (at richain = 1) as a function of chain length. This value always remains positive indicating stability against phase separation. 82 3.4 Summary The compression of a single end-grafted chain by a disc of finite radius was investigated using numerical SCF theory. Chain segment distributions for uncompressed chains (mushrooms) are consistent with results predicted by de Gennes [1980]. However, in contrast to more recent work using a model based on scaling theory [Subramanian et al. 1995; 1996] we find no first- order escape transition with hysteresis when a chain is compressed. Based on our model, the chain appears to be squeezed out from under the disc gradually as the disc approaches the surface. The energy required to compress the mushroom increases monotonically with compression and becomes independent of chain length following chain escape. 83 4 Brushes and Brush-Particle Interactions 4.1 Outline of Calculations Figure 4.1 Schematic diagram showing a solute molecule interacting with polymer chains end-grafted to a surface. The centerline represents the axis of a cylindrical coordinate system, and the rings on the surface represent spatial increments in the radial direction. Because of applications for end-grafted polymer surfaces in the field of biomaterials and bioseparations, the interactions between solute molecules and the end-grafted polymer chains becomes important to understand. Several molecular-based models have been developed for analyzing brush distributions as a function of conditions. Few studies however, have focused on a small impenetrable particle (a protein molecule perhaps) interacting with, or moving into, polymer chains that are end-grafted to a surface. In this work we have studied the interaction between a cylindrical particle and a brush using our cylindrical SCF model (Figure 4.1). The density distribution of an undisturbed brush has been calculated using different interaction parameters, and as a function of surface density and chain length. The brush density distribution around a particle has also been calculated for varying system conditions and displacement of polymer segments by the particle has been monitored. Finally, the excess energy to move the particle into the brush has been evaluated. 84 4.2 Brush Density Distributions A thorough discussion of brush density distributions has been presented in Chapter 1. Here, we provide a few examples of results from our model, which are in excellent agreement with experiment [Cosgrove, 1990a; 1990b] and provide a significantly more detailed and accurate picture than previous analytical theories. Figure 4.2 shows density profiles for a 50 segment brush (b) at 5% grafting density under various solvent (o) conditions. The polymer- solvent interaction parameter Xbo has been varied from -0.5 to 1.0. The polymer-surface (s) and solvent-surface interaction is athermal (i.e. %os - %ps = 0). Layer(z) Figure 4.2 Brush density distributions for various polymer-solvent interaction parameters. N= 50 segments, o = 5%, j o s = %\>* = 0. A more favourable brush-solvent interaction (i.e. more negative) leads to a more extended polymer layer. The brush profile for the case of j b o = 10 is almost completely 85 collapsed, forming a very dense surface layer. In a worse than theta solvent (#,0 > 0.5 as N - » QO) the polymer layer collapses since free polymer in solution would prefer to phase separate under these conditions. The brush rapidly expands away from the grafting surface as Xbo decreases from poor {xbo = 1) to athermal (xbo = 0) solvent conditions. However, further improvement of the solvent {xbo < 0) does not promote as significant an enhancement in chain stretching. In a miscible solvent, the degree of chain stretching is determined by a balance between the loss in chain entropy when extended to a linear conformation, and the gain in polymer-segment stepping paths when adjacent lattice sites become occupied by solvent. The energy required to stretch a chain into a more linear state scales with brush height squared, while the energy gain accompanying an increase in polymer-segment free volume scales approximately linearly with brush height. Thus, the stretching penalty becomes prohibitive at large brush heights and results in a very weak dependence of brush height on %00 when %bo < 0. The effect of polymer-surface interactions on the brush density profile is shown in Figure 4.3. For chains of 50 segments at 5% surface density, the polymer-surface interaction determines whether or not a depletion zone is seen next to the grafting surface. For a strong polymer-surface attraction, the depletion zone disappears. Although the overall brush profile is not changed much at high enough surface densities, a strongly adsorbing surface can lead to a profile that decreases exponentially away from the surface at low surface densities (pancake). 86 Layer(z) Figure 4.3 Brush density profiles for various polymer-surface interactions are shown. N = 50 segments, a= 5%, %os = %bo = 0. Brush density distributions are shown in Figure 4.4 for four different chain lengths at 10% surface density with %b0 = 0.5. As the chain length is increased, the brush extends further away from the surface. Relatively little increase is observed in the maximum density, especially at the longer chain lengths. Both results are in agreement with analytical predictions [Milner et al. 1988a; de Gennes, 1976]. The brush height increases approximately linearly with chain length. In contrast to analytical model results [Milner et al., 1988a] the segment density in the outer edge of the brush is predicted to decrease more slowly, forming an extended tail region. Due to this effect, it is common to discuss the brush height in terms of the root mean squared layer thickness, given by - jz2<b(z)dz j|<D(z)dz. It is also worth noting that the depletion and tail regions occupy a smaller fraction of the overall distribution as chain length increases. 87 O(Z) Layer (z) Figure 4.4 Segment density profiles for brushes of different lengths. Surface density is 10% and the polymer-solvent interaction parameter b̂o is 0.5. All other interactions are athermal. Layer (z) Figure 4.5 Brush density profiles for various surface densities. Chain length is 100 segments and %bo = 0.5. 88 Figure 4.5 shows brush density distributions for 100 segment chains at four different surface grafting densities. Increasing surface density causes the maximum density of the brush to increase, as scaling models predict [Milner, 1991]. Also the increased excluded volume interactions lead to a more stretched brush and hence the height scales (less than linearly) with surface density. 4.3 Brush-Particle Interactions 4.3.1 2D Brush Density Distributions Introduction of a particle changes the conformations of chains in the brush (Figure 4.6). The distribution of chains of N = 50 segments at a grafting density of cr =0.1 is shown in Figure 4.6 for an interacting particle of radius Rp = 3 and length Lp = 3 (referred to as a 3 by 3 particle) at various distances from the surface. In order to try and simulate physiological conditions, we have set %bo , %po , and Xbp to 0.4, 0.4 and 0, respectively, for most of our calculations. In addition, the calculations involve non-surface-mobile chains, unless otherwise specified. All other Xij are set to zero. These values were chosen based on the interaction parameters calculated from the LALLS data of Haynes et al. [1993] discussed in Chapter 2. In this way, the results should provide at least a semi-quantitative picture of globular protein interactions with grafted PEG brushes, which are finding increased application as biomaterials and anti-fouling surfaces. It is worth noting, however, that the energy curves due to brush interaction with the particle are not very sensitive to the brush-solvent or particle-solvent interaction parameters, although they are sensitive to a favorable interaction between the brush and particle. 89 Figure 4.6 Polymer segment density distributions in the cylindrical lattice are indexed in both the z and r directions. The chains of length N=50 with a grafting density of cr=0.1 are interacting with a particle of dimensions Rp=3 and Lp=3 (3 by 3). Interaction parameters are set as Xbo=0.4 and Xpo=0-4 (unspecified interaction parameters are always zero). The density distribution is shown for a particle position (a) outside the brush, (b) at layer zp = 10, (c) zp = 7, and (d) zp = 5. A brush density profile is shown in Figure 4.6 (a) at conditions where the particle is positioned outside the brush. The distribution is completely uniform in the radial direction, and shows a typical brush profile in the z-direction at near theta solvent conditions. The z and r coordinates on the plot index the mean-field radial of the cylinder, whose average brush segment volume fraction is shown on the vertical axis. In Figure 4.6 (b), (c) and (d), brush distributions are shown when the particle is sitting at layers zp = 10, zp = 7, and zp= 5, respectively. For lattice sites occupied by the particle, the volume fraction of solvent and brush is essentially zero. 90 There is a slight decrease in the brush density near the edges of the particle due to the entropic penalty of chains stepping next to an impenetrable object. As the particle begins to penetrate the brush, the segment density fills in quickly behind the particle, indicating chain splaying. To reach this extended chain state there is a concomitant depletion of the segment density between the particle and the surface (in front of the particle). When the particle is positioned in layer 5 (Figure 4.6 (d)), the segment density profile also increases slightly at the sides of the particle. This leads to a physical picture of grafted chains spread around the particle such that they radially partition away from the centerline, thereby lowering the brush density in the center of the particle's path. Segment depletion between the particle and grafting surface does not occur at all conditions, however. Figure 4.7 shows the brush density distribution for four different chain lengths (at constant surface density), with a 3 by 3 particle sitting at layer 5. For a chain length of 15 segments, Figure 4.7 (a) shows that the segment density is slightly higher between the particle and the surface, indicating brush compression without significant chain splaying. Segment depletion is first observed at a chain length of N =25 and becomes increasingly more pronounced with increasing chain length. Thus, the configuration of the grafted chains in the presence of a penetrating particle is a strong function of the chain length relative to the particle radius. This relationship has not been reported previously since the modeling methods based on the Alexander-de Gennes brush [Alexander, 1977; de Gennes, 1976] and analytical SCF theory [Milner et al., 1988a] do not allow splaying to occur. 91 Figure 4.7 The 2-dimensional segment density distributions for four different chain lengths are shown around a 3 by 3 particle placed with the leading edge at zp = 5. Interaction parameters are specified as Xbo=0.4 and Xpo=0.4. Chain lengths are (a) N = 15, (b) N = 25, (c) N = 40 and (d) N = 60 segments. Changing the surface density of the grafted chains also affects the density distribution of the brush. Figure 4.8 shows a graft of N = 50 segment chains interacting with a 5 by 5 particle placed at layer 7. The chains are grafted at densities of <j = 0.05, cr = 0.1, cr = 0.2, and cr = 0.3 in Figure 4.8 (a), (b), (c), and (d), respectively. Increasing the surface density results in an extension of the brush along with an increase in the maximum in its density distribution. At high surface densities, the brush profile begins to approach a step-profile in accordance with the Alexander-de Gennes brush model. Unlike variations in chain length (see Figure 4.7), increasing 92 the grafting density does not in general lead to a transition from complete compression of the brush to a regime where chain splaying makes a dominant contribution to the perturbed density distribution. Instead, it changes the relative amount of segment compression or splaying as seen in the segment density distribution. As the surface density increases, changes in segment density between the particle and the surface due to particle approach become less significant with respect to the rest of the density profile. Figure 4.8 The brush distribution of N=50 segment chains is shown around a 5 by 5 particle at four different grafting densities. Interaction parameters are Xbo=0.4 and Xpo=0.4. Grafting densities are (a) a=0.05, (b) G=0.1, (C) CT=0.2 and (d) a=0.3. 93 4.3.2 Brush-Particle Interaction Energy The excess energy A'n< required to move different sized particles into a brush is shown in Figure 4.9. For a brush made up of 50-segment chains at a graft density of cr =0.1, A'"' increases with increasing particle size, particularly aCdeeper penetration depths. Under the chosen conditions, formation of contacts between brush and particle segments is enthalpically favorable since both segments have a net unfavorable interaction with the solvent (xto >0). Nevertheless, the calculated A'"' values are always positive and increase sharply as the particle penetrates deeper into the brush. Thus, A'"' is dominated by entropic effects, which are strongly repulsive due to the chain conformations that are eliminated by the presence of the particle. Layer Position (zp) Figure 4.9 Brush-particle interaction energy curves are shown for chains of N=50 with a graft density of a=0.1 interacting with various sized particles. Interaction parameters are specified as Xbo-0.4 and Xpo=0.4. Particle size has been varied from Rp=l, Lp=l up to Rp=7, Lp=7. The dashed line shows the inflection point of each curve that shows one. 94 The effect of chain length on brush-particle interaction energies is shown in Figure 4.10 for a 3 by 3 particle. As the chains are made longer the onset of an energetic repulsion occurs farther from the surface. This distance appears to scale linearly with N, as we might expect since h ~ N [Milner, 1991]. For longer chains, Amt becomes nearly independent of chain length at particle positions close to the surface. This might be explained by the fact that increased chain length extends the brush, and allows chains to fill in the space behind the particle more easily. Once a particle has significantly penetrated a relatively long brush, the conformational entropy lost by the chains due to particle approach is at least partly offset by the gain of configurations of the chain segments, which step behind the particle. Figure 4.10 Interaction energy curves are shown for a particle of Rp=3, Lp=3 interacting with different length polymer chains of constant grafting density a - 0.1. Interaction parameters are Xbo=0.4 and Xpo=0.4. The dashed line shows the inflection point of each curve. 95 The dotted lines in Figure 4.9 and Figure 4.10 intersect the energy curves at a point of inflection. Each inflection point represents a maximum in the force F required to move the particles into the brush, since F - -dAml Idzp. An inflection point is not observed for all combinations of 7Y and Rp. Instead, there is a critical chain length Ncmcai above which an inflection point is seen for a given particle size. Figure 4.11 shows Ncritical as a function of Rp for a constant grafting density of 10%. The critical chain length for this system increases linearly with Rp. A linear dependence of the critical chain length on Rp corresponds to a linear increase in the distance that chain segments must reach in order to occupy space behind the particle. The minimum chain length for an inflection point in the limit of an infinitesimal particle is ca. 5 segments. 80 60 I 40- 20- 0 Ncrttical = 4.9 + 9.5 Rp 3 R. 5 Figure 4.11 The critical chain length Ncriticai for which an inflection point is observed is shown as a function of Rp for the system with a = 0.1, the previous interaction parameters, and particles with a cylindrical geometry where Lp = RD. 96 The change in brush-segment density profile which leads to the observed inflection point was analyzed by calculating A'"' curves for penetration of a 3 by 3 particle under two limiting cases: (i) pure compression of the brush by the particle (which is equivalent to the model of Subramanian et al. [1996]), and (ii) penetration into a brush in which all chains grafted directly beneath the particle have been removed (so that no chain splaying can occur). Placing the reflecting boundary of the lattice in radial Rp + 1, where Rp is the particle radius, allows the energy of pure chain compression to be calculated. Under this limiting condition, A'"' increases monotonically with particle approach and no inflection point is observed. Thus, the inflection point is related to chain segment escape from beneath the approaching particle. In principle, the inflection point could result from the lower energy of splayed chains relative to compressed chains and/or the increase in configurational entropy of chains which extend beyond the particle volume and can thereby step into the solvent-rich space behind the particle. Results from limiting case (ii), where no chain splaying can occur, show an inflection point, suggesting that the inflection is due to the invasion of extended chain segments into the space behind the particle. Indeed, in this case A'"' begins to decrease (indicating a slight attraction) as the particle is moved very close (zp < 3) to the grafting surface. Segment density profiles under this condition indicate that the chains adjacent to the side walls of the particle reach a lower energy state by stretching slightly in the z-direction to allow a larger number of segments to sample the volume directly behind the particle. A similar result, in which there is an inflection point and a decrease in A'"' when the particle is close to the surface, is obtained when model calculations are carried out for a brush of surface-mobile chains (section 4.3.3). It is also worth noting that the model of Subramanian et al., which is limited by the severe assumption that 97 chain splaying can not occur, does not predict the presence of an inflection point. As a result, Subramanian et al. predict the existence of a maximum in compression force only for a surface covered by multiple mobile mushrooms. 1.5-, 1.0 2 o 0.5- 0.0 (0 4 -c E cu 2 o m Q_ CO Q 0 c (D E <D 0 10 15 20 25 Layer Position (z ) Figure 4.12 F and the rate of segment displacement from the total volume the particle moves through are shown as a function of zp for a 3 by 3 particle moving into a brush with N = 50 and a = 0.1. Our results indicate that an inflection point can appear when chain segments are able to fill in the space behind the particle. As it penetrates the brush, the particle reduces the conformational entropy of the brush, increasing A'"'. When the chains are able to gain a significant amount of conformational entropy by filling in behind the particle, the rate of increase in A'"' with respect to particle position gets smaller, giving rise to the inflection point. It is therefore useful to consider the segment density in the volume defined by the particle and the volume directly above and below it. As the particle is moved toward the surface, the number of 98 chain segments in the particle's path volume is always decreasing. As segments begin to fill in behind the particle, however, the rate of segment displacement changes. Figure 4.12 shows the rate of decrease in the total number of segments (in front of and behind the particle) as a function of particle position; -dA,nt/dzp is also shown for the same system. In systems where N is above Ncritical, the inflection point in A'"' always coincides with the point where the rate of displacement of segments by the particle is decreasing the fastest with respect to zp. That is, where the rate of accumulation of segments behind the particle is a maximum. In the MD simulations by Murat and Grest [1996], the segment density between an AFM tip and the surface was found to remain approximately constant due to rearrangement of the polymer chains around the tip. The radius of their AFM cylinder was large compared to the segment length (7 to 16 times) for simulations with 100 segment chains. Under such conditions, we also see little or no depletion of segments if particle size is made large. They also found a constantly increasing repulsive force under compression. Our results are consistent with the fact that they used a long cylinder, which extends outside the brush to model the AFM tip. As a result, there was no empty space behind the tip that would allow chains segments to fill in behind it. The dependence of the inflection point position Zinflection and the magnitude of the force maximum Fmax on particle size, chain length and surface density were also investigated for specific cases. For a brush with N = 50 and a= 0.1 The inflection point was found to move toward the surface in roughly linear proportion with increasing Rp (Figure 4.9). For the same system, Fmax was found to increase in proportion to R2*, as shown in Figure 4.13, which suggests that it scales closely with particle volume. For a 3 by 3 particle and 10% grafting density (Figure 4.10), Zinflecti0n moves away from the surface almost linearly with increasing N 99 (^inflection ~ N0'8), which corresponds closely to the scaling of brush height. Furthermore, Fmax was found to decrease with increasing chain lengths such that Fmax ~ N ~1'2. Longer brush chains can more quickly adopt configurations that step behind a particle. This combined with the fact that a particle of fixed size will take away a smaller fraction of the conformations sampled by larger chains, results in an earlier inflection point in A ' n t and lower Fmax. max 2.5 2.0 1.5- 1.0- 0.5- 0.0 F ~ R max p - i 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 R Figure 4.13 Fmax is shown as a function of particle size for a brush with N= 50 and a 10% graft density. The maximum force appears to scale closely with particle volume. The dependence of Amt, Zj„flection and Fmax on surface density was also investigated for a brush with N- 50 and 3 by 3 particle and is shown in Figure 4.14. The energetic repulsion was found to increase significantly with grafting density, due primarily to the increased overall density in the brush distribution. The onset of repulsion also appears to start at a distance from the surface that increases linearly with grafting density. This agrees with the observation that 100 brush height scales linearly with a in a very poor solvent (x = 0-4) [Fleer et al., 1993]. The values of Zp and Fmax were both found to increase approximately linearly with the grafting density (Z,^ec,;o„ ~ a °'8, Fmax. ~ cr). — CJ=0.3 • a =0.25 a =0.2 -—CT=0.15 — a=0.1 ---a=0.05 a =0.02 15 20 25 30 Layer Position (zp) Figure 4.14 Interaction energy curves are shown for a particle of Rp=3, Lp=3 interacting with 50 segment polymer chains of varied grafting density. Interaction parameters are Xbo=0.4 and Xpo = 0.4. The dashed line shows the inflection point of each curve. 4.3.3 Surface Mobile Brushes Grafting can be simulated such that the chain segments attached at layer 1 redistribute themselves radially when a particle is placed very close to the surface, simulating what might occur if the graft were surface mobile. The chain ends can also be fixed such that each grafted segment is confined to a specific radial at z=\. In the latter case, the brush is fixed to the grafting surface by individually grafting chains to each radial. The number of chains used to calculate the 101 normalization constant for each separate radial must then be calculated using the number of lattice sites in the radial and assuming a constant grafting density. The fixed brush therefore consists of Rmax components, for each of which Equation (2.39) has a different starting point. Finally the total brush density distribution is calculated by summing the volume fractions of the separate components. Figure 4.15 Brush distribution around a 3 by 3 particle is shown for surface non- mobile and mobile brushes. Brushes consist of N=40 segment chains at 10% surface density and interaction parameters used were Xbo=0.4 and Xpo-0.4. The distributions for non-mobile chains around a particle in (a) layer zp - 4 and (b) zp = 2 are shown; distributions of mobile chains around a particle in (c) zp - 4 and (d) zp = 2 are also shown. A comparison of the segment density distribution for mobile and non-mobile chains around a particle is shown in Figure 4.15. The brush distributions were calculated for N = 40 segment chains at 10% surface density around a 3 by 3 particle. Figure 4.15 (a) and (b) show 102 non-mobile chains with the particle positioned in layer zp - 4 and zp = 2, respectively. The non- mobile chains can not decrease their density in layer z = 1 below the grafting density of 10% (0 = 0.1). As a result, the chain segment density at the maximum can be seen to increase around the edges of the particle. Figure 4.15 (c) and (d) show the density distribution of the surface mobile brush around a particle in layer zp = 4 and 2, respectively. The brush density for the mobile chains at the maximum is uniform in the radial distribution indicating that the chains are able to redistribute themselves within the lattice. When the particle is positioned at zp = 2, the density between the particle and the surface in Figure 4.15 (d) is significantly below 0 = 0.1, indicating that chains have been displaced from the volume directly beneath the particle. Figure 4.16 Interaction energy curves are shown for a 3 by 3 particle interacting with 40 segment polymer chains that are surface mobile (dashed) and non- mobile (solid). Interaction parameters are %bo=0.4 and Xpo=0.4. 103 The brush-particle interaction energies A'"' for the mobile and non-mobile brushes are compared in Figure 4.16 (for the same conditions as Figure 4.15). At close approach to the surface (zp < 6) the mobile brush density directly beneath the particle becomes lower than that of the non-mobile brush, thus providing a weaker entropic repulsion. At a particle position directly next to the surface (zp= 1) the interaction energy actually decreases for the mobile brush due to an increase in the configurational entropy of the chains. This suggests that if the particle were able to penetrate all the way into the brush it would prefer to remain there. 4.3.4 Brush-Particle Attraction Joen et al. [1991a; 1991b] have predicted a minimum in the interaction energy (negative A'"') between a protein and a polyethylene oxide graft when a strong hydrophobic interaction free energy term is included in their model. On closer approach of the particle to the grafting surface, they predict a strong steric repulsion that leads to a monotonically increasing interaction free energy. In order to compare our model predictions with those of Joen et al. we have set the brush-solvent Xbo, particle-solvent Xpo, and the brush-particle Xbp interaction parameters to 0.5, 0.5 and -0.5 respectively (all other Xy = 0). These interaction parameters will be assumed for the remainder of this discussion. A strongly negative brush-particle interaction parameter is needed in order to see a minimum in A m t with our model. The magnitude of XbP is almost unrealistic for an interaction that is not electrostatic in nature, suggesting that perhaps Joen et al. might have over-estimated the hydrophobic interaction energy. Experiments have shown that short PEG chains can increase protein adsorption slightly [Merril, 1992; Kishida et al., 1992] suggesting a net weak attraction between polymer and protein. Interaction parameters calculated from data of Haynes et al. 104 [1993] show that the polymer-protein interaction is weak, however, with a maximum attractive value of XbP = -0.1. The interaction energy for various sized particles with a 50-segment brush at 10% grafting density is shown in Figure 4.17. As the particles begin to move into the brush, favorable brush-particle energetic contacts dominate the interaction energy, leading to negative values of A'"'. With further penetration, a net energetic repulsion occurs once the loss of chain entropy is greater than the favorable contact energy. As the particle size increases, the position of the onset of energetic repulsion Z o m e t occurs sooner (higher zp). As shown in Figure 4.17, the increase in Zomet with particle size gets progressively smaller as the particle size becomes larger. This suggests that Z o n s e t will reach a constant value, equal to that for compression by an infinite flat plate. 5 10 15 20 Layer Position (zp) Figure 4.17 Brush-particle interaction energy curves are shown for five different sized particles interacting with chains of length N=50 and grafting density ci=0.1. Interaction parameters are Xbo^O.S, Xpo=0.5 and XbP= -0.5. 105 The smallest particle (1 by 1) occupies ca. 3 lattice sites and because of its small size does not experience a significant entropic repulsion by the chains. As a result A'"' is negative for the 1 by 1 particle beyond layer 2 due to the energetic attraction. For the particles with Rp > 1, as Rp increases, the magnitude of the attractive minimum \A^\ increases significantly such that |^min|~^/> 7 S (ignoring the 1 by 1 particle). This favorable interaction energy scales closely with particle surface area, since it is based on the number of contacts that brush and solvent segments make with the surface of the particle. The minimum in Ami indicates that the solute particles can preferentially adsorb to the surface of a brush under appropriate conditions. CT =0.05 CT =0.1 CT =0.15 CT =0.2 0 5 10 15 20 25 Layer Position (z ) Figure 4.18 Brush-particle interaction energy curves are shown for a 5 by 5 particle penetrating chains of length N=50 at four different grafting densities. Interaction parameters are %\,o=0.5, %vo=0.5 and XbP= -0.5. 106 The dependence of A'" on cr is shown in Figure 4.18. The brush consists of 50 segment chains that are interacting with a 5 by 5 particle. As the grafting density is increased, Zomet increases with increasing <j(Zonset ~ cr04) due to the extension of the brush with increased graft density. The model predicts that A^n decreases slightly with increasing cr and is shifted away from the grafting surface. The gradual increase in the attractive minimum may result from more brush-particle contacts at higher surface density. The entropic repulsion appears to scale more strongly with increasing a than the energetic attraction, and as a result the attractive regions of the Aml curves become narrower with increasing cr. These two effects agree with the observation that a brush profile becomes more step-like as surface density is increased. However, the dependence of Al^m on a is quite weak and the scaling therefore difficult to determine. We have also investigated the effect of N on the interaction energy for this set of interaction parameters as shown in Figure 4.19. Zomet shifts away from the surface approximately linearly with increasing N, again corresponding to the extension of the brush from the surface. The minimum in A'"' becomes slightly shallower with increasing N, although the dependence appears very weak. However, longer chains have a broader attractive minimum than short chains. This may be because long chains are able to form multiple contacts with the particle over a wider range of particle positions. If the criterion for protection of a surface is the distance at which a repulsion energy begins, Zomet, then for a given brush length and surface density, we predict that larger particles should be repelled better. Also, because h and Zomet increase with increasing N and cr, maximum repulsion will be obtained with long chains at the highest surface density possible. 107 0 5 10 15 20 25 30 35 Figure 4.19 Brush-particle interaction energy curves are shown for a 3 by 3 particle penetrating brushes of varied chain length, each at 10% grafting density. Interaction parameters are xt>o=0.5, Xpo=0.5 and X b P = -0.5. The magnitude and width of the attractive minimum may also be important considerations if it is important to minimize any possible interaction between a protein and the brush. For example, studies have shown that although grafted polymer will strongly inhibit platelet adhesion, transient contacts can still lead to platelet activation and decreased platelet counts [Llanos and Sefton, 1993]. If the criteria for choosing the optimum surface coverage are maximizing Zomet while at the same time minimizing the concentration of particles in the region of attractive interaction with the brush, predicting the optimum conditions becomes more difficult. In the ideal solution limit, the Boltzmann factor Kp(zp) = exp(-y4mt (zp)/kT), gives the concentration of solute particles Cp(zp) at position zp with respect to the bulk solution as Cp(zp) = 108 Cp" * Kp(zp). The region of negative interaction energy will therefore have a higher concentration of particles than in the bulk solution, while the repulsive region has a much lower concentration. In order to determine optimum brush conditions when considering both the magnitude and width of the minimum, it is useful to calculate the average concentration of solute particles in the brush. For simplicity, we assume that the brush is an infinite plane in the x and y- directions (perpendicular to the grafting surface). The average particle concentration in the brush is then given by c - » = c j * v f _ ( 4 1 ) K where h is the brush height arbitrarily chosen by us as the point at which the volume fraction of an undisturbed brush falls below 0.1%. Multiplying the average concentration by the height of the brush gives a value that is proportional to the number of particles that are interacting with the brush. Figure 4.20 shows Kp as a function of particle position zp for different sized particles interacting with a brush of 50 segments. At positions zp where Kp is greater than unity (negative A1"'), the concentration of solute particles will be greater than that in bulk solution. Based on Equation (4.1), the average particle concentration in the brush C^s is equal to the area under each curve (up to the brush height h) divided by h. The smallest particle (1 by 1) will partition strongly into the brush due to the attractive energetic interaction and a lack of entropic repulsion. For the larger particles, however, the entropic repulsion is much more significant and Kp is seen 109 to fall below unity upon any appreciable penetration into the brush. The attractive region at the outer edge of the brush increases with particle size as seen by the increase in Kp. We predict that the total number of particles interacting with the brush decreases with increasing a due to the narrower attractive region at higher a (Figure 4.18). Therefore, increasing surface density for a given chain length should always improve resistance to solute interaction with a surface, even in cases where brush-particle attraction occurs. Figure 4.20 Layer partition coefficient Kp(zp) for various particle sizes interacting with 50 segment chains at 10% surface density. Interaction parameters are Xbo=0.5, xPo=0.5 and Xb P= -0.5. When chain length is increased, the width of the attractive minimum is increased causing an increase in partitioning of particles to the brush surface region (Figure 4.19). This suggests that very long flexible chains may actually be less effective at reducing adsorption, especially since experimentally, lower surface densities are normally obtained with higher degrees of 110 polymerization. Under conditions of weak attraction, the optimum properties of end-grafted polymer for protection against adsorption are therefore a combination of low N and high <r that gives an acceptable value of ZonSet- It is important to note that we find a different result for optimum protein repulsion than Joen et al. [1991a; 1991b]. Their steric repulsion term, based on an energy balance between a chain stretching term and an osmotic energy term, does not allow the brush to splay around the protein. We have shown that splaying is an important phenomenon in the case where a small particle penetrates a brush. Their model assumes a constant density for the brush, which increases as the protein moves toward the surface. As a result, they overestimate the density of polymer between the particle and the surface and thereby predict very large repulsive energies. By adding a steric repulsion term to a hydrophobic interaction energy term, they have assumed that the terms are thermodynamically additive. This is almost certainly a poor assumption, since entropy plays a significant role in the "steric" repulsion by the polymer chains as well as in the hydrophobic interaction energy. Other work which considers the mechanism of protein repulsion by polymers such as PEG predicts that it may be due to the hydrophilic nature of the polymer, or the strong hydration shell required to solvate the polymer [Gblander et al., 1992]. Our calculations suggest that while those may contribute, the dominant effect in protein repulsion is the high degree of flexibility in the polymer chains, which in thermodynamic terms is seen as high conformational entropy. Penetration of a particle into the chains decreases the chain entropy, giving rise to a sharp increase in the brush-particle interaction energy. This is supported by the results of Schroen et al. [1995] who found that an extended brush configuration of adsorbed Pluronic copolymer was 111 significantly better at reducing protein adsorption than a collapsed layer, even though the surface was made more hydrophilic in both cases. 112 NOMENCLATURE Symbols a = activity (mol/L, g/L) = experimental constant used to calculate radius of gyration (nm) A = Helmholtz energy Au, A i U = osmotic second (mL-mol/g2) and third virial coefficient (mL2-mol/g3) A m o n = maximum monolayer capacity of grafted polymer (g/m2) B = proportionality constant for calculation of analytical brush density profile c = weight concentration of solute (g/mL) = conformation of a polymer chain in the lattice model C, = normalization constant for volume fraction calculation in numerical SCF model Cp(zp) = concentration of the particle at position zp within the brush Cp = average concentration of solute within the brush D = distance between two grafting surfaces upon brush compression ft = probability that any site adjacent to a previously occupied site is already filled during placement of chain i+1 in a lattice / = vector set of nonlinear functions used to solve numerical SCF model F = sum of squares of nonlinear functions used to measure the extent of solution convergence G(z, r) = free segment weighting factor; weight for a free segment to be in lattice position (z,r) with respect to the bulk solution. G(z, r, s) = chain end distribution function; weight for segment s of the chain to be in lattice position (z,r) h = brush height J = Jacobian matrix used for solution of nonlinear set of equations k = Boltzmann constant (1.3 8x10"23 J/K) 113 Kp = Boltzmann factor giving the ratio of solute concentration inside the brush to that in bulk solution (partition coeffient) /, IL = polymer repeat unit length; lattice unit length for one chain segment (nm) If. = number of lattice sites available for placement of the first chain segment of the first chain in the cylindrical lattice L(r) = number of lattice sites in one layer of radial r of the cylindrical lattice m = mass of one polymer molecule Mm = mass of one monomer unit (kDa) M, M„, Mw = polymer molecular weight; number and weight averaged molecular weight rij = number of molecules of component i present in the solution or lattice ni = number of chains of type /" adopting a conformation defined by c n(r) = total number of polymer segments in radial r summed over all lattice layers z N = number of segments in a lattice chain = degree of polymerization of a polymer molecule NB, NLB = number of bonds in a polymer molecule; number of segment bonds in a lattice chain p, PL = polymer chain persistence (a measure of stiffness); lattice chain persistence based on lattice geometry P° = pure component vapour pressure (kPa) Q = canonical partition function (fixed N, V, and T) r = radial position index in cylindrical lattice geometry r. (z, r) = number of segments of a chain of type /' adopting conformation c that are in lattice radial (z,r) R = root mean squared end-to-end distance of a polymer molecule in solution Rg = radius of gyration; root mean squared distance of segments from the center of mass of a polymer molecule in solution Rmax = maximum radial dimension in cylindrical lattice 114 Rrms = root mean squared distance of chain segments from the cylindrical axis s = segment ranking number S(r) = surface area of the outer edge of radial r (lattice units squared) S, Sm, Sid = entropy; entropy of mixing; ideal entropy of mixing (J/K) T = absolute temperature (K) u = potential of mean force u' = hard core potential (space filling potential) umt = interaction potential u = vector of iteration variables consisting of potential of mean force values U = internal energy V] = molar volume of solvent (L/mol) Xi = mole fraction of component /' z = layer position index in lattice model = maximum layer position in cylindrical lattice Greek Letters <E> = volume fraction T = experimental surface concentration (g/m2) IT = osmotic pressure Q = multiplicity term; total possible number of ways of arranging the system S = grand canonical partition function (fixed p, V, and T) a(z,r) = Lagrange multiplier that ensures lattice sites in position (z,r) are filled completely, but without segment overlap % = Flory interaction parameter 115 X, X = lattice step probability; product of step probabilities for a chain adopting conformation c ju = chemical potential v = proportionality constant in analytical SCF theory vl+i = number of sites available for placing chain i+l into a lattice already containing /' chains v = partial specific volume of a dissolved solute (mL/g) 9 = surface coverage in equivalent monolayers (or total number of grafted chain segments per surface lattice site) p = density (g/mL) a = lattice model surface grafting density or surface fraction (chains/surface site) < H c = number of ways of placing a single volumeless chain of type / with conformation c into the lattice a>ij = energy associated with contact between one segment of /" and one segment of/ Acoij = interchange energy; difference in contact energy between one unlike contact [12] and half the sum of two different like contacts [11], [22] 116 REFERENCES 1. Alexander, S. (1977) Adsorption of Chain Molecules with a Polar Head: A Scaling Description. J. Phys. (Paris), 38, 977-987. 2. Andrade, J.D. and V. Hlady (1986) Protein Adsorption and Materials Biocompatibility: A Tutorial Review and Suggested Hypotheses. Advances in Polymer Science (K. Dusek, ed.) Springer-Verlag: Berlin, 1-63. 3. Ansarifar, M.A. and P.F. Luckham (1988) Measurement of the Interaction Force Profiles Between Block Copolymers of Poly(2-vinylpyridine)/Poly(t-butylstyrene) in a Good Solvent. Polymer, 29, 329-335. 4. Auroy, P., L. Auvray and L. Leger (1991a) Structures of End-Grafted Polymer Layers: A Small-Angle Neutron Scattering Study. Macromolecules. 24,2523-2528. 5. Auroy, P., L. Auvray and L. Leger (1991b) Characterization of the Brush Regime for Grafted Polymer Layers at the Solid-Liquid Interface. Physical Review Letters, 66(6), 719- 722. 6. Bamford, C.H. and K.G. Al-Lamee (1992) Chemical Methods for Improving the Haemocompatibility of Synthetic Polymers. Clinical Materials, 10, 243-261. 7. Baranowski, R. and M.D. Whitmore (1997) Theory of the Structure of Adsorbed Block Copolymers: Detailed Comparison with Experiment. J. Chem. Phys., 103(6), 2343-2353. 8. Baskir, J.N., T.A. Hatton and U.W. Suter (1989a) Affinity Partitioning in Two-Phase Aqueous Polymer Systems: A Simple Model for the Distribution of the Polymer-Ligand Tail Segments near the Surface of a Particle. J. Phys. Chem., 93(2), 969-976. 9. Baskir, J.N., T.A. Hatton and U.W. Suter (1989b) Thermodynamics of the Partitioning of Biomaterials in Two-Phase Aqueous Polymer Systems: Comparison of Lattice Model to Experimental Data. J. Phys. Chem., 93, 2111-2122. 10. Baskir, J.N., T.A. Hatton and U.W. Suter (1987) Thermodynamics of the Separation of Biomaterials in Two-Phase Aqueous Polymer Systems: Effect of the Phase-Forming Polymers. Macromolecules. 20, 1300-1311. 11. Bayer, E. and W. Rapp (1992) Polystyrene-Immobilized PEG Chains: Dynamics and Application in Peptide Synthesis. Immunology, and Chromatography. Poly(Ethylene Glycol) Chemistry: Biotechnical and Biomedical Applications (J.M. Harris, ed.) Plenum Press: New York, 325-345. 117 12. Bergstrom, K., K. Holmberg, A. Safranj, A.S. Hoffman, M.J. Edgell, A. Kozlowski, B.A. Hovanes and JM. Harris (1992) Reduction of Fibrinogen Adsorption on PEG-coated Polystyrene Surfaces. J. Biomedical Materials Research, 26, 779-790. 13. Brandrup, J. and E.H. Immergut (1989) Polymer Handbook, Wiley: New York. 14. Brent, R.P. (1973) Algorithms for Minimization Without Derivatives, Prentice-Hall Inc.: New York. 15. Brooks, D.E. and V.J. Miiller (1996) Size Exclusion Phases and Repulsive Protein-Polymer Interaction/Recognition. J. Molecular Recognition, 9, 697-700. 16. Broyden, C.G. (1965) A Class of Methods for Solving Nonlinear Simultaneous Equations. Mathematical Computing, 19, 577-593. 17. Chakrabarti, A., P. Nelson and R.J. Toral (1994) Interpenetrations in Polymer Brushes. J. Chem. Phys., 100, 748-749. 18. Chakrabarti, A. and R. Toral (1990) Density Profile of Terminally Anchored Polymer Chains: A Monte Carlo Study. Macromolecules. 23, 2016-2021. 19. Claesson, P.M. et al. (1996) Techniques for Measuring Surface Forces. Advances in Colloid and Interface Sci. ,67, 119-183. 20. Clayfield, E.J. and E.C. Lumb (1966) A Theoretical Approach for Polymeric Dispersant Action I. Calculation of Entropic Repulsion Exerted by Random Polymer Chains Terminally Adsorbed on Plane Surfaces and Spherical Particles. J. Colloid Interface Sci., 22, 269-284. 21. Cohen Stuart, M.A. (1980) Ph.D. Thesis, Wageningen Agricultural University. 22. Cosgrove, T., T.G. Heath, J.S. Phipps and RM. Richardson (1991) Neutron Reflectivity Studies of Polymers Adsorbed on Mica from Solution. Macromolecules, 24, 94-98. 23. Cosgrove, T. (1990a) Volume-Fraction Profiles of Adsorbed Polymers. J. Chem. Soc. Faraday Trans., 86(9), 1323-1332. 24. Cosgrove, T. and K. Ryan (1990b) NMR and Neutron-Scattering Studies on Polyfethylene oxide) Terminally Attached at the Polystyrene/Water Interface. Langmuir, 6, 136-142. 25. Cosgrove, T., T. Heath, B. van Lent, F. Leermakers and J. Scheutjens (1987a) Configuration of Terminally Attached Chains at the Solid/Solvent Interface: Self- Consistent Field Theory and a Monte Carlo Model. Macromolecules, 20, 1692-1696. 118 26. Cosgrove, T., T. Heath, G. Ryan and B. van Lent (1987b) The Conformation of Adsorbed Poly(ethylene oxide) at the Polystyrene/Water Interface. Polymer Communications, 28, 64- 65. 27. Davis, S.S. (1997) Biomedical Applications of Nanotechnology - Implications for Drug Targeting and Gene Therapy. Trends in Biotechnology, 15, 217-224. 28. de Gennes, P.-G. (1985) CR Acad. Sci. Paris, Ser. II, 300, 839. 29. de Gennes, P.-G. (1980) Conformations of Polymers Attached to an Interface. Macromolecules, 13, 1069-1075. 30. de Gennes, P.-G. (1976) Scaling Theory of Polymer Adsorption. J. Phys. (Paris), 37, 1443- 1452. 31. Dennis, J.E. and R.B. Schnabel (1983) Numerical Methods for Unconstrained Optimization and Nonlinear Equations. Prentice-Hall, Inc.: Englewood Cliffs, N.J. 32. Desai, N.P. and J.A. Hubbell (1991) Biological Responses to Polyethylene Oxide Modified Polyethylene Terephthalate Surfaces. J. Biomedical Materials Research, 25, 829-843. 33. DiMarzio, E.A. and R.J. Rubin (1971) Adsorption of a Chain Polymer between Two Plates. J. Chem. Phys., 55(9), 4318-4336. 34. Doroszkowski, A. and R. Lambourne (1971) Measurement of the Strength of Steric Barriers in Non-Aqueous Polymer Dispersions. J. Polymer Sci.: Part C, 34, 253-264. 35. Du, H., P. Chandaroy and S.W. Hui (1997) Grafted Polyethylene glycol) on Lipid Surfaces Inhibits Protein Adsorption and Cell Adhesion. Biochimica et Biophysica Acta, 1326, 236-248. 36. Dunn, S.E., A. Brindley, S.S. Davis, M.C. Davies and L. Ilium (1994) Polvstvrene- Polv(Ethvlene Glycol) (PS-PEG2000) Particles as Model Systems for Site Specific Drug Delivery. 2. The Effect of PEG Surface Density on the in Vitro Cell Interaction and in Vivo Biodistribution. Pharmaceutical Research. 11(7), 1016-1022. 37. Elbert, D.L. and J.A. Hubbell (1996) Surface Treatments of Polymers for Biocompatibility. Annu. Rev. Mater. Sci., 26, 365-394. 38. Evers, O.A., J.M.H.M. Scheutjens and G.J. Fleer (1991) Statistical Thermodynamics of Block Copolymer Adsorption. 3. Interaction between Adsorbed Layers of Block Copolymers. Macromolecules, 24, 5558-5566. 119 39. Evers, O.A., J.M.H.M. Scheutjens and G.J. Fleer (1990) Statistical Thermodynamics of Block Copolymer Adsorption. 1: Formulation of the Model and Results for the Adsorbed Layer Structure. Macromolecules. 23, 5221-5233. 40. Field, J.B., C. Toprakcioglu, R.C. Ball, HB. Stanley, L. Dai, W. Barford, J. Penfold, G. Smith and W. Hamilton (1992) Determination of End-Adsorbed Polymer Density Profiles by Neutron Reflectometry. Macromolecules. 25, 434-439. 41. Fleer, G.J., M.A. Cohen Stuart, J.M.H.M. Scheutjens, T. Cosgrove, B. Vincent (1993) Polymers at Interfaces, Chapman & Hall: London. 42. Flory, P. (1969) Statistical Mechanics of Chain Molecules, John Wiley & Sons, Inc.: New York. 43. Flory, P. (1953) Principles of Polymer Chemistry, Cornell University Press: Ithaca, N.Y. 44. Flory, P. (1942) Thermodynamics of High Polymer Solutions. J. Chem. Phys., 10, 51-61. 45. Golander, C.G., J.N. Herron, K. Lim, P. Claesson, P. Stenius and J.D. Andrade (1992) Properties of Immobilized PEG Films and the Interaction with Proteins: Experiments and Modeling. Poly(Ethylene Glycol) Chemistry: Biotechnical and Biomedical Applications (J.M. Harris, ed.) Plenum Press: New York, 221-245. 46. Gombotz, W.R., W. Guanghui, T.A. Horbett and A.S. Hoffman (1992) Protein Adsorption to and Elution from Polyether Surfaces. PolyfEthylene Glycol) Chemistry: Biotechnical and Biomedical Applications (J.M. Harris, ed.) Plenum Press: New York, 247-261. 47. Gombotz, W.R., W. Guanghui, T.A. Horbett and A.S. Hoffman (1991) Protein Adsorption to Poly(ethylene oxide) Surfaces. J. Biomedical Materials Research, 25, 1547-1562. 48. Grest, G.S. and M. Murat (1993) Structure of Grafted Polymeric Brushes in Solvents of Varying Quality: A Molecular Dynamics Study. Macromolecules, 26, 3108-3117. 49. Guffond, M.C., D.R.M. Williams and EM. Sevick (1997) End-Tethered Polymer Chains under AFM Tips: Compression and Escape in Theta Solvents. Langmuir, 13(21), 5691- 5696. 50. Hadziioannou, G., S. Patel, S. Granick and M. Tirrell (1986) Forces Between Surfaces of Block Copolymers Adsorbed on Mica. J. American Chemical Society, 108, 2869-2876. 51. Haynes, CA. and W. Norde (1994) Globular Proteins at Solid/Liquid Interfaces. Colloids and Surfaces B: Biointerfaces, 2, 517-566. 120 52. Haynes, C.A., F.J. Benitez, H.W. Blanch, J.M. Prausnitz (1993) Application of Integral- Equation Theory to Aqueous Two-Phase Partitioning Systems. AIChE Journal, 39, 1539- 1557. 53. Haynes, C.A., R.A. Beynon, R.S. King, H.W. Blanch and J.M. Prausnitz (1989) Thermodynamic properties of Aqueous Polymer Solutions: Poly(ethylene glycoD/Dextran, J. Phys. Chem., 93, 5612-5617. 54. Hesselink, F.T. (1971) On the Theory of the Stabilization of Dispersions by Adsorbed Macromolecules. I. Statistics of the Change of Some Configurational Properties of Adsorbed Macromolecules on the Approach of an Impenetrable Interface. J. Phys. Chem., 75(1), 65-71. 55. Huang, K. and A.C. Balazs (1993) A Two-Dimensional Self-Consistent-Field Model for Grafted Chains: Determining the Properties of Grafted Homopolymers in Poor Solvents. Macromolecules, 26, 4736-4738. 56. Huggins, M.L. (1942) Some Properties of Solutions of Long Chain Compounds. J. Phys. Chem., 46, 151-158. 57. Jeon, S.I., J.H. Lee, J.D. Andrade and P.-G. de Gennes (1991a) Protein-Surface Interactions in the Presence of Polyethylene Oxide I. Simplified Theory. J. Colloid Interface Sci., 142(1), 149-158. 58. Jeon, S.I and J.D. Andrade (1991b) Protein-Surface Interactions in the Presence of Polyethylene Oxide I. Effect of Protein Size. J. Colloid Interface Sci., 142(1), 159-166. 59. Johansson, H.O., G. Karlstrom, F. Tjerneld and CA. Haynes (1997) Driving Forces for Phase Separation and Partitioning in Aqueous Two-Phase Systems, accepted January, 1997 to J. Chromatography B, in press. 60. Kent, M.S., L.T. Lee, B.J. Factor, F. Rondelez and G.S. Smith (1995) Tethered Chains in Good Solvent Conditions: An Experimenal Study Involving Langmuir Diblock Copolymer Monolayers. J. Chem. Phys., 103(6), 2320-2342. 61. Kent, M.S., L.T. Lee, B. Farnoux and R. Rondelez (1992) Characterization of Diblock Copolymer Monolayers at the Liquid-Air Interface by Neutron Reflectivity and Surface Tension Measurements. Macromolecules, 25, 6240-6247. 62. Kenworthy, A.K., K. Hristova, D. Needham and T.J. Mcintosh (1995) Range and Magnitude of the Steric Pressure Between Bilayers Containing Phospholipids with Covalently Attached Poly(ethylene glycol). Biophysical Journal, 68, 1921-1936. 121 63. King, R.S., H.W. Blanch and J.M. Prausnitz (1988) Molecular Thermodynamics of Aqueous Two-Phase Systems for Bio separations, AIChE Journal, 34(10), 1585-1594. 64. Kishida, A., K. Mishima, E. Corretge, H. Konishi and Y. Ikada (1992) Interactions of Poly(ethylene glycoD-grafted Cellulose Membranes with Proteins and Platelets, Biomaterials, 13(2), 113-118. 65. Kuhn, W. (1934) KolloidZ, 68, 2. 66. Lai, P.Y. and K. Binder (1992) Structure and Dynamics of Polymer Brushes Near the 8 Point: A Monte Carlo Simulation. J. Chem. Phys., 97(1), 586-595. 67. Lai, P.Y. and K. Binder (1991) Structure and Dynamics of Grafted Polymer Layers: A Monte Carlo Simulation. J. Chem. Phys., 95(12), 9288-9299. 68. Lee, L.T. and M.S. Kent (1997) Direct Measurement of the Penetration of Free Chains into a Tethered Chain Laver. Physical Review Letters, 79(15), 2899-2902. 69. Lee, L.T., B.J. Factor, M.S. Kent and F. Rondelez (1994) J. Chem: Soc. Faraday Discuss., 98, 130. 70. Leermakers, F.A.M., J.M.H.M. Scheutjens and J. Lyklema (1990) Statistical Thermodynamics of Association Colloids. IV. Inhomogeneous Membrane Systems, Biochimica et Biophysica Acta, 1024, 139-151. 71. Leermakers, F.A.M. and J.M.H.M. Scheutjens (1988) Statistical Thermodynamics of Association Colloids. I. Lipid Bilaver Membranes. J. Chem. Phys., 89(5), 3264-3274. 72. Lim, K. and J.N. Herron (1992) Molecular Simulation of Protein-PEG Interaction. Poly(Ethylene Glycol) Chemistry: Biotechnical and Biomedical Applications (J.M. Harris, ed.) Plenum Press: New York, 29-56. 73. Llanos, G.R. and M.V. Sefton (1993) Immobilization of Poly(ethylene glycol) onto a Poly(vinyl alcohol) Hydrogel: 2. Evaluation of Thrombogenicity. J. Biomedical Materials Research, 27, 1383-1391. 74. Luckham, P.F. (1996) Recent Advances in Polymers at Surfaces: The Steric Effect, Current Opinion in Colloid and Interface Science, 1, 39-47. 75. Malmsten, M. and J.M. Van Alstine (1996) Adsorption of Poly(Ethylene Glycol) Amphiphiles to Form Coatings Which Inhibit Protein Adsorption, J. Colloid Interface Sci., Ill, 502-512. 122 76. Martin, J.I. and Z.G. Wang (1995) Polymer Brushes: Scaling. Compression Forces. Interbrush Penetration, and Solvent Size Effects. J. Phys. Chem., 99, 2833-2844. 77. McMillan, Jr., W.G. and J.E. Mayer (1945) The Statistical Thermodynamics of Multicomponent Systems. J. Chem. Phys., 13, 276. 78. Meier, D.J. (1967) Theory of Polymeric Dispersants. Statistics of Constrained Polymer Chains. J. Phys. Chem., 71(6), 1861-1868. 79. Merrill, E.W. (1992) PolyfEthylene Oxidê and Blood Contact. A Chronicle of One Laboratory. Poly(Ethylene Glycol) Chemistry: Biotechnical and Biomedical Applications (J.M. Harris, ed.) Plenum Press: New York, 199-220. 80. Milner, S.T. (1991) Polymer Brushes. Science, 251, 905-914. 81. Milner, S.T., T.A. Witten and M.E. Cates (1988a) Theory of the Grafted Polymer Brush. Macromolecules, 21, 2610-2619. 82. Milner, S.T., T.A. Witten and M.E. Cates (1988b) A Parabolic Density Profile for Grafted Polymers. Europhysics Letters, 5(5), 413-418. 83. Murat, M. and G.S. Grest (1996) Molecular Dynamics Simulations of the Force Between a Polymer Brush and an AFM Tip. Macromolecules, 29, 8282-8284. 84. Murat, M. and G.S. Grest (1991) Interaction Between Grafted Polymeric Brushes. Computer Simulation of Polymers (R.J. Roe, ed.) Prentice-Hall Inc.: Englewood Cliffs, N. J., 141-153. 85. Murat, M. and G.S. Grest (1989) Structure of a Grafted Polymer Brush: A Molecular Dynamics Simulation. Macromolecules, 22, 4054-4059. 86. Napper, D.H. (1983) Polymeric Stabilization of Colloidal Dispersions, Academic Press: New York. 87. Patel, S., M. Tirrell and G. Hadziioannou (1988) A Simple Model for Forces Between Surfaces Bearing Grafted Polymers Applied to Data on Adsorbed Block Copolymers. Colloids and Surfaces, 31, 157-179. 88. Prausnitz, J.M., R.N. Lichtenthaler and E.G. de Azevedo (1986) Molecular Thermodynamics of Fluid-Phase Equiliibria, 2nd ed., Prentice-Hall Inc.: Englewood Cliffs, N.J. 123 89. Privalov, P.L. (1979) Stability of Proteins: Small Globular Proteins. Adv. Prot. Chem., 33, 167. 90. Rathbone, S.J., C.A. Haynes, H.W. Blanch and J.M. Prausnitz (1990) Thermodynamic Properties of Dilute Aqueous Polymer Solutions from Low-Angle Laser-Light-Scattering Measurements. Macromolecules. 23, 3944-3947. 91. Ruckenstein, E. and B.Q. Li (1997) Steric Interactions between Two Grafted Polymer Brushes. J. Chem. Phys., 107(3), 932-942. 92. Scales, L.E. (1985) Introduction fo Non-Linear Optimization, Springer-Verlag: New York. 93. Scheutjens, J.M.H.M. and G.J. Fleer (1985) Interaction between Two Adsorbed Polymer Layers. Macromolecules. 18, 1882-1900. 94. Scheutjens, J.M.H.M. and G.J. Fleer (1980) Statistical Theory of the Adsorption of Interacting Chain Molecules. 2. Train. Loop and Tail Size Distribution. J. Phys. Chem., 84, 178-190. 95. Scheutjens, J.M.H.M. and G.J. Fleer (1979) Statistical Theory of the Adsorption of Interacting Chain Molecules. 1. Partition Function. Segment Density Distribution, and Adsorption Isotherms. J. Phys. Chem., 83, 1619-1635. 96. Schroen, C.G.P.H., M.A. Cohen Stuart, K. van der Voort Maarschalk, A. van der Padt and K. van't Reit (1995) Influence of Preadsorbed Block Copolymers on Protein Adsorption: Surface Properties. Layer Thickness, and Surface Coverage. Langmuir, 11, 3068-3074. 97. Semenov, A.N. (1985) Contribution to the Theory of Microphase Layering in Block Copolymer Melts. JEW, 61, 733. 98. Silberberg, A. (1968) Adsorption of Flexible Macromolecules IV. Effect of Solvent-Solute Interactions. Solute Concentration and Molecular Weight. J. Chem. Phys., 46, 2835-2851. 99. Subramanian, G., D.RM. Williams and P.A. Pincus (1996) Interaction between Finite- Sized Particles and End Grafted Polymers. Macromolecules, 29, 4045-4050. 100. Subramanian, G., D.R.M. Williams and P.A. Pincus (1995) Escape Transitions and Force Laws for Compressed Polymer Mushrooms. Europhysics Letters, 29(4), 285-290. 101. Taunton, H.J., C. Toprakcioglu, L.J. Fetters and J. Klein (1990) Interactions between Surfaces Bearing End-Adsorbed Chains in a Good Solvent. Macromolecules. 23, 571-580. 124 102. Taunton, H.J., C. Toprakcioglu, L.J. Fetters and J. Klein (1988) Forces between Surfaces Bearing Terminally Anchored Polymer Chains in Good Solvents. Nature, 332, 712-714. 103. Tyn, M.T. and T.W. Gusek (1990) Prediction of Diffusion Coefficients of Proteins. Biotechnology andBioengineering, 35, 327-338. 104. Uchida, E. and Y. Ikada (1997) Topography of Polymer Chains Grafted on a Polymer Surface Underwater. Macromolecules, 30, 5464-5469. 105. Van Alstine, J.M. and M. Malmsten (1997) Poly(ethylene glycol) Amphiphiles: Surface Behavior of Biotechnical Significance. Langmuir, 13, 4044-4053. 106. van Lent, B., R. Israels, J.M.H.M. Scheutjens and G.J. Fleer (1990) Interaction between Hairy Surfaces and the Effect of Free Polymer. J. Colloid Interface Sci., 137, 380-394. 107. Wijmans, CM. and B.J. Factor (1996) Self-Consistent Field Calculations of Tethered Chains in the Presence of Mobile Homopolymer: Comparison with Experiment. Macromolecules, 29, 4406-4411. 108. Wijmans, CM., E.B. Zhulina and G.J. Fleer (1994a) Effect of Free Polymer on the Structure of a Polymer Brush and Interaction between Two Polymer Brushes. Macromolecules, 27, 3238-3248. 109. Wijmans, CM., F.A.M. Leermakers and G.J. Fleer (1994b) Pair Potentials Between Polymer-Coated Mesoscopic Particles. Langmuir, 10, 4514-4516. 110. Wijmans, CM., F.A.M. Leermakers and G.J. Fleer (1994c) Chain Stiffness and Bond Correlations in Polymer Brushes. J. Chem. Phys., 101(9), 8214-8223. 111. Wijmans, CM., J.M.H.M. Scheutjens and E.B. Zhulina (1992) Self-Consistent Field Theories for Polymer Brushes. Lattice Calculations and an Asymptotic Analytical Description. Macromolecules, 25, 2657-2665. 112. Williams, D.R.M. and F.C MacKintosh (1995) Polymer Mushrooms Compressed Under Curved Surfaces. J. Phys. II France, 5, 1407-1417. 113. Woodle, M.C and D.D. Lasic (1992) Sterically Stabilized Liposomes. Biochimica et Biophysica Acta, 1113, 171-199. 114. Yang, Z. and H. Yu (1997) Preserving a Globular Protein Shape on Glass Slides: A Self- Assembled Monolayer Approach. Advanced Materials. 9(5), 426-429. 125 115. Yeung, C, A.C. Balazs and D. Jasnow (1993) Lateral Instabilities in a Grafted Layer in a Poor Solvent. Macromolecules, 26, 1914-1921. 116. Zhulina, E.B., V.A. Pryamitsyn and O.V. Borisov (1989) Structure and Conformational Transitions in Grafted Polymer Chain Layers. A New Theory. Polymer Science U.S.S.R., 31(1), 205-216. 117. Zubay, G.L. (1993) Biochemistry, Wm. C. Brown Publishers: Dubuque, Iowa. 126 Appendix A Fortran Code for Cylindrical SCF Model 127 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * PROGRAM SCF_CYLINDER ********************************************************************* C C This program c a l c u l a t e s the volume f r a c t i o n d i s t r i b u t i o n of a C g r a f t e d polymer, solvent and a g r a f t e d p a r t i c l e i n a C c y l i n d r i c a l geometry according to the Scheutjens and F l e e r C model. C C Reference: Scheutjens, F l e e r , Cosgrove, Vincent, Cohen St u a r t , C "Polymers at I n t e r f a c e s " , Cambridge: Chapman and C H a l l , 1993 C C By Brad S t e e l s C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C p **** L i s t o f V a r i a b l e s : MAIN **** c ISEG_B LENGTH OF BRUSH MOLECULES c ISEG_P LENGTH OF PARTICLE MOLECULES c BRUSH NUMBER OF BRUSH MOLECULES c PARTICL NUMBER OF PARTICLE MOLECULES c I LAY NUMBER OF LAYERS IN LATTICE c IRAD NUMBER OF RADIALS IN LATTICE c NSEGS NUMBER OF DIFFERENT SEGMENT TYPES IN SYSTEM c M TOTAL NUMBER OF VARIABLES (ILAY*IRAD*NSEGS) c CHI BO BRUSH-SOLVENT INTERACTION PARAM. c CHIBS BRUSH-SURFACE INTERACTION PARAM. c CHIBP BRUSH-PARTICLE INTERACTION PARAM. c CHIPO PARTICLE-SOLVENT INTERACTION PARAM. c CHIPS PARTICLE-SURFACE INTERACTION PARAM. c U(i) MEAN POTENTIAL (ITERATION VARIABLE) c F ( i ) FUNCTION VALUE c ENTR(i) CONFORMATIONAL ENTROPY c ENGY(i) INTERACTION ENERGY c TOTENG(i) HELMHOLTZ ENERGY CHANGE W.R.T. UNMIXED COMPONENTS c IP0S1, IPOS2, IPOS c p START, FINISH AND CURRENT POSITION OF PARTICLE I — c **** Subroutine NEWTON **** c DFDU(i,j) JACOBIAN MATRIX c H ( i , j ) INVERSE JACOBIAN c FNEW(i) NEW FUNCTION VALUES c UNEW(i) NEW POTENTIAL VALUES c DF(i) CHANGE IN FUNCTION VALUES FROM LAST ITERATION c DU(i) CHANGE IN POTENTIAL VALUES FROM LAST ITERATION c p ALPHA FRACTIONAL STEP TAKEN IN SEARCH DIRECTION c C **** Subroutine CYLINDER **** c L NUMBER OF LATTICE SITES IN A GIVEN RADIAL c EQUAL TO PI*(2R-1), R BEING INTEGER POSITION c L0,L1 LATTICE STEP PROBABILITIES IN Z DIRECTION 128 C LR,LR0,LR1 : RADIAL STEP PROBABILITIES (DEPEND ON R) C PHIP ( z,r), PHIB ( z,r), PHIO(z,r) C : VOLUME FRACTION OF PARTICLE, BRUSH AND SOLVENT C GEND(z,r,s) : CHAIN END WEIGHTING FACTORS C GTERM(z,r,S): CHAIN END WEIGHTING FACTORS -- GRAFTED CHAIN C UOHARD(z,r), UBHARD(z, r ) , UPHARD(z,r) C : HARD CORE POTENTIAL FOR SOLVENT, BRUSH AND C PARTICLE C ************************************************ IMPLICIT REAL*8(A-H,0-Z) PARAMETER(MD = 2 0 00,MAXZ = 5 0,MAXR=2 0) DIMENSION PHIB(MAXZ,MAXR), PHIO(MAXZ,MAXR), PHIP(MAXZ,MAXR), 1 U(MD), ENTR(MAXZ), ENGY(MAXZ), TOTENG(MAXZ) COMMON/PARAMS/CHIBO,CHIBP,CHIBS,CHIPO,CHIPS,CHIOS,BRUSH, 1 PARTICL,NSEGS COMMON/POSITION/IPOS1,IPOS2,IPOS C OPEN (UNIT=14,FILE='OUT.DAT') OPEN (UNIT=1,FILE='LOG' ) C CALL INPUTS(ILAY,IRAD,ISEG_B,ISEG_P,M) C MAXIT = 1000 C START_TIME = F_TIME( ) C C Loop to move p a r t i c l e p o s i t i o n . C DO I = IPOS2, IPOS1, -1 IPOS = I C CALL INITIALIZE(M,ILAY,IRAD,U) C CALL NEWTON(M,U,ILAY,IRAD,ISEG_B,ISEG_P,MAXIT, 1 PHIB,PHIP,PHIO) C CALL ENERGY(U,ILAY,IRAD,ISEG_B,ISEG_P,PHIO,PHIB,PHIP, 1 ENTR(I),ENGY(I),TOTENG(I)) C FINAL = (F_TIME( ) - START_TIME) I_MIN = INT(FINAL/60) I_SEC = INT(FINAL - I_MIN*60) C CALL OUTPUT(ILAY,IRAD,ISEG_B,ISEG_P,U,PHIB,PHIP,PHIO, 1 I_MIN,I_SEC,M,ENTR,ENGY,TOTENG) C END DO C CLOSE(UNIT=14) CLOSE(UNIT=1) C END C 129 c c * * c c * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * SUBROUTINE INPUTS(ILAY,IRAD,ISEG_B,ISEG_P,M) Parameters necessary f o r execution are read from the input f i l e and p r i n t e d to the screen and output f i l e . ******************************************************************* IMPLICIT REAL*8(A-H,0-Z) COMMON/PARAMS/CHIBO,CHIBP,CHIBS,CHIPO,CHIPS,CHIOS,BRUSH, 1 PARTICL, NSEGS COMMON/POSITION/IPOS1,IPOS2,IPOS OPEN ( UNIT = 10, FILE = ' i n l . d a t ' ) C C READ (10, *' NSEGS READ (10, *' ISEG_B READ(10, *' ISEG_P READ(10, *' ILAY READ(10, *' IRAD READ(10, *' IPOS1 READ (10, *' IPOS2 READ(10, *' CHIBO READ(10, *' CHIBS READ(10, *' CHIBP READ(10, *' CHIPO READ (10, *' CHIPS CHIOS = 0 DO READ(10, * SURDEN READ(10, *' PARTICL !START POSITION OF PARTICLE CENTER !END POSITION OF PARTICLE CENTER !PERCENT SURFACE DENSITY CLOSE(UNIT=10) M = ILAY * IRAD * NSEGS BRUSH = SURDEN / 1.D2 * 3.14159D0 * DBLE(IRAD)**2.DO WRITE(6,'(2X, A39,14 )') 1 'VALUE FOR GRAFTED CHAIN LENGTH WRITE(6,'(2X,A39,14)') 1 'VALUE FOR PARTICLE CHAIN LENGTH WRITE(6,'(2X,A39,I4)') 1 'VALUE FOR CYLINDRICAL LATTICE LENGTH WRITE(6,'(2X,A39,I4)') 1 'VALUE FOR CYLINDER RADIUS WRITE(6,'(2X,A39,14)') 1 'START POSITION FOR PARTICLE WRITE(6,'(2X,A39,14)') 1 'END POSITION FOR PARTICLE WRITE(6,'(2X,A39,G12.4)') 1 'BRUSH-SOLVENT INTERACTION PARAMETER WRITE(6, 1(2X,A39,G12.4)') 1 'BRUSH-SURFACE INTERACTION PARAMETER WRITE(6, ' (2X,A3 9,G12.4) ') 1 'BRUSH-PARTICLE INTERACTION PARAMETER WRITE(6,'(2X,A39,G12.4)') 1 'PARTICLE-SOLVENT INTERACTION PARAMETER WRITE(6,'(2X,A39,G12.4)') , ISEG_B , ISEG_P , ILAY , IRAD , IPOS1 , IPOS2 , CHIBO , CHIBS , CHIBP , CHIPO 130 1 'PARTICLE-SURFACE INTERACTION PARAMETER WRITE(6, '(2X,A39,G12.4)') 1 'SURF DENS. OF BRUSH MOLECULES WRITE(6, ' (2X,A39(G12.4) ') 1 'NUMBER OF BRUSH MOLECULES WRITE(6,'(2X,A39,G12.4)') 1 'NUMBER OF PARTICLE MOLECULES WRITE(6, ' (2X,A39,I4) ' ) 1 'NUMBER OF ITERATION VARIABLES WRITE(14,'(2X, A39,I4 )') 1 'VALUE FOR GRAFTED CHAIN LENGTH WRITE(14, ' (2X,A39,14) ' ) 1 'VALUE FOR PARTICLE CHAIN LENGTH WRITE(14,'(2X,A39,14)') 1 'VALUE FOR CYLINDRICAL LATTICE LENGTH WRITE(14,'(2X,A3 9,14)') 1 'VALUE FOR CYLINDER RADIUS WRITE(14,'(2X,A39,I4)') 1 'START POSITION FOR PARTICLE WRITE(14,'(2X,A39,14)') 1 'END POSITION FOR PARTICLE WRITE(14, ' (2X,A39,G12.4) ' ) 1 'BRUSH-SOLVENT INTERACTION PARAMETER WRITE(14, 1(2X,A3 9,G12.4)') 1 'BRUSH-SURFACE INTERACTION PARAMETER WRITE(14,'(2X,A39,G12.4)') 1 'BRUSH-PARTICLE INTERACTION PARAMETER WRITE(14,'(2X,A39,G12.4)') 1 'PARTICLE-SOLVENT INTERACTION PARAMETER WRITE(14,'(2X,A39,G12.4)') 1 'PARTICLE-SURFACE INTERACTION PARAMETER WRITE(14,'(2X,A39,G12.4)') 1 'SURF DENS. OF BRUSH MOLECULES WRITE(14,'(2X,A39,G12.4)') 1 'NUMBER OF BRUSH MOLECULES IN LATTICE WRITE(14, ' (2X,A39,G12.4) ' ) 1 'NUMBER OF PARTICLE MOLECULES WRITE(14, ' (2X,A39,I4) ' ) 1 'NUMBER OF ITERATION VARIABLES CHIPS SURDEN BRUSH PARTICL M ISEG_B ISEG_P ILAY IRAD IPOS1 IPOS2 CHIBO CHIBS CHIBP CHIPO CHIPS SURDEN BRUSH PARTICL M RETURN END * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * SUBROUTINE NEWTON(M,U,ILAY,IRAD,ISEG_B,ISEG_P,MAXIT, 1 PHIB,PHIP,PHIO) ********************************************************************* IMPLICIT REAL*8(A-H,0-Z) PARAMETER(MD=2000,MAXZ=50,MAXR=20) EXTERNAL FUNLIN DIMENSION U(MD), F(MD), DFDU(MD,MD), FNEW(MD), UNEW(MD), DU(MD), 1 DF(MD), H(MD,MD), PHIB(MAXZ,MAXR), PHIO(MAXZ,MAXR), 2 PHIP(MAXZ,MAXR) 131 c TOL = l.D-7 ITERS = 0 C C I n i t i a l i z e parameters f o r the l i n e search r o u t i n e . C XI = l.D-3 XF = 7.D-1 XTOL = l.D-5 FTOL = l.D-6 IMAX = 5 C C Carry out an i n i t i a l f u n c t i o n e v a l u a t i o n . C CALL CYLINDER(U,ILAY,IRAD,ISEG_B,ISEG_P,F,PHIB,PHIP,PHIO) SUMSQR = SSQ(F,M) C WRITE(6,*) WRITE(6,*) ' CALCULATING JACOBIAN ' C CALL NUM_DRIV(F,U,ILAY,IRAD,ISEG_B,ISEG_P,M,PHIB,PHIP, 1 PHIO,DFDU) C WRITE(6,*) WRITE(6,*) ' CALCULATING INVERSE OF JACOBIAN ' C CALL MATINV(DFDU,M,H,IERROR) C WRITE(6,*) WRITE(6,*) ' SEARCHING FOR SOLUTION OF NONLINEAR SYSTEM ' WRITE(6,*) WRITE(6,*) ' ITERS ',' SSQ (old) (new) ',' ALPHA' WRITE(1,*) ' ITERS ',' SSQ (old) (new) ',' ALPHA' C DO WHILE ( SQRT(SUMSQR) .GT. TOL ) ITERS = ITERS + 1 IF (ITERS .GT. MAXIT) THEN RETURN END IF C C Obtain the search d i r e c t i o n DU from H and F. C DO 20 I = 1, M SUM = 0.D0 DO 10 J = 1, M SUM = SUM - H(I,J) * F(J) 10 CONTINUE DU(I) = SUM 2 0 CONTINUE C DO 25 I = 1, M FNEW(I) = F(I) UNEW(I) = U(I) 25 CONTINUE 132 c C Carry out the l i n e a r search i n the d i r e c t i o n DU. C CALL LINSRCH(FUNLIN,XI,XF,XTOL,FTOL,IMAX,FNEW,UNEW,PHIB,PHIP, 1 PHIO,ILAY,IRAD,ISEG_B,ISEG_P,M,DU,SSQNEW,X,IERR) C WRITE(6,30) ITERS, SUMSQR, SSQNEW, X WRITE(1,30) ITERS, SUMSQR, SSQNEW, X 30 FORMAT(I5,3(G15.5)) C C Update the f u n c t i o n and v a r i a b l e values us i n g X (alpha). C DO 4 0 I = 1, M DF(I) = FNEW(I) - F(I) F(I) = FNEW(I) U(I) = UNEW(I) DU(I) = X * DU(I) 4 0 CONTINUE SUMSQR = SSQNEW C C Update the inverse Jacobian using the new f u n c t i o n values C CALL APROX_DRIV(DF,DU,M,H) C END DO RETURN END C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * DOUBLE PRECISION FUNCTION FUNLIN(X,U,ILAY,IRAD,ISEG_B,ISEG_P, 1 M,DU,F,PHIB,PHIP,PHIO) C Function c a l l e d by LINSRCH to minimize SSQ w.r.t. X. * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * IMPLICIT REAL*8(A-H,O-Z) PARAMETER(MD=2000,MAXZ=50,MAXR=20) DIMENSION DU(MD), U(MD), F(MD), UU(MD), PHIB(MAXZ,MAXR), 1 PHIO(MAXZ,MAXR), PHIP(MAXZ,MAXR) COMMON/LINEAR/UU C C Update the temporary v a r i a b l e UU and t r y a f u n c t i o n e v a l u a t i o n . C DO 10 I = 1, M UU(I) = U(I) + X * DU(I) 10 CONTINUE CALL CYLINDER(UU,ILAY,IRAD,ISEG_B,ISEG_P,F,PHIB,PHIP,PHIO) FUNLIN = SSQ(F,M) C RETURN END C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * SUBROUTINE LINSRCH(FUNC,XI,XF,XTOL,FTOL,IMAX,F,UP,PHIB,PHIP, 1 PHIO,ILAY,IRAD,ISEG_B,ISEG_P,NVAR,DU,Y,X,IERR) C C a r r i e s out a u n i v a r i a t e minimization on the sum of squares 133 C with respect to the step f r a c t i o n , X, i n the search d i r e c t i o n . * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * IMPLICIT REAL*8(A-H,0-Z) REAL*8 M PARAMETER(MD=2000,MAXZ=50,MAXR=20) DIMENSION F(MD), UP(MD), DU(MD), UU(MD), PHIB(MAXZ,MAXR), 1 PHIP(MAXZ,MAXR), PHIO(MAXZ,MAXR) COMMON/LINEAR/UU EXTERNAL FUNC C C I n i t i a l i z e the search i n t e r v a l on A and B. C X has the l e a s t value of F, W the next lowest and V the C previous value of W. U i s the current e v a l u a t i o n p o i n t . C A = XI B = XF C = (3.DO - DSQRT(5.D0)) / 2.DO V = A + C * (B-A) W = V X = V E = 0.D0 FV = FUNC(V,UP,ILAY,IRAD,ISEG_B,ISEG_P,NVAR,DU,F,PHIB, 1 PHIP,PHIO) FW = FV FX = FV IFEVAL = 0 IERR = 0 M = 0.D0 T2 = 0.D0 C C Check the stopping c r i t e r i a f o r a l o c a l minimum. C DO WHILE ( DABS(X-M) .GT. (T2-(.5D0*(B-A))) ) C IFEVAL = IFEVAL + 1 C C Check f o r maximum number of i n t e r p o l a t i o n s . C IF (IFEVAL .GT. IMAX) THEN IERR = 1 Y = FX DO 10 I = 1, NVAR UP(I) = UU(I) 10 CONTINUE RETURN END IF C M = .5D0 * (A+B) TOL = XTOL * DABS(X) + FTOL ! R e l a t i v e and absolute c r i t e r i a T2 = 2 * TOL P = 0.D0 Q = P R = P IF ( DABS(E) .GT. TOL ) THEN ! P a r a b o l i c I n t e r p o l a t i o n 134 R = (X-W) * (FX-FV) Q = (X-V) * (FX-FW) P = (X-V) * Q - (X-W) * R Q = 2 * (Q-R) IF ( Q .GT. O.DO ) THEN P = -P ELSE Q = -Q END IF R = E E = D END IF IF ( DABS(P) .LT. DABS(.5D0*Q*R) .AND. P .GT. (Q*(A-X)) 1 .AND. P .LT. (Q*(B-X)) ) THEN D = P / Q U = X + D ! Take a quadratic step C C Cannot evaluate F too clos e to A or B. C IF ( (U-A) .LT. T2 .OR. (B-U) .LT. T2 ) THEN D = TOL IF ( X .GE. M ) D = D * (-1.D0) END IF ELSE IF ( X .LT. M) THEN ! Golden S e c t i o n step taken E = B - X ELSE E = A - X END IF D = C * E C C Cannot evaluate F too clos e to X. C IF ( DABS(D) .GE. TOL ) THEN U = X + D ELSEIF ( D .GT. O.DO ) THEN U = X + TOL ELSE U = X - TOL END IF END IF FU = FUNC(U,UP,ILAY,IRAD,ISEG_B,ISEG_P,NVAR,DU,F,PHIB,PHIP,PHIO) C C Update the values of a, b, v, w, and x. C IF ( FU .LE. FX ) THEN IF ( U .LT. X ) THEN B = X FB = FX ELSE A = X FA = FX END IF V = W 135 FV = FW W = X FW = FX X = U FX = FU ELSE IF ( U .LT. X ) THEN A = U FA = FU ELSE B = U FB = FU END IF IF ( FU .LE. FW .OR. W .EQ. X ) THEN V = W FV = FW W = U FW = FU ELSEIF ( FU .LE. FV .OR. V .EQ. X .OR. V .EQ. W ) THEN V = U FV = FU END IF END IF END DO Y = FX DO 100 1 = 1 , NVAR UP(I) = UU(I) 100 CONTINUE C RETURN END C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * SUBROUTINE APROX_DRIV(DF,DU,M,H) C Update the inverse Jacobian using Broyden's update. *************************************************************** IMPLICIT REAL*8(A-H,O-Z) PARAMETER (MD=2 000) DIMENSION H(MD,MD), DU(M), DF(M), HDF(MD), DUH(MD) DO 2 0 I = 1, M SUM = 0.D0 DO 10 J = 1, M SUM = SUM + H(I,J) * DF(J) 10 CONTINUE HDF(I) = SUM 20 CONTINUE DENOM = 0.D0 DO 4 0 I = 1, M SUM = 0.D0 DO 30 J = 1, M SUM = SUM + DU(J) * H(J,I) 3 0 CONTINUE DUH(I) = SUM DENOM = DENOM + SUM * DF(I) 136 4 0 CONTINUE DO 8 0 I = 1, M DO 70 J = 1, M ADD = (HDF(I) - DU(I)) * DUH(J) / DENOM H(I,J) = H(I,J) - ADD 70 CONTINUE 8 0 CONTINUE RETURN END C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * SUBROUTINE NUM_DRIV(F,U,ILAY,IRAD,ISEG_B,ISEG_P,M,PHIB,PHIP, 1 PHIO,DFDU) C D e r i v a t i v e s are approximated by f i n i t e d i f f e r e n c e s . A r e l a t i v e C step s i z e , DELU, of about SQRT(machine p r e c i s i o n ) * U i s used and C i f U i s zero, then a step s i z e of SQRT(machine p r e c i s i o n ) i s C used. In double p r e c i s i o n , machine p r e c i s i o n = 2.2e-16 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * IMPLICIT REAL*8(A-H,0-Z) PARAMETER(MD=2000,MAXZ=50,MAXR=20) DIMENSION U(MD), DFDU(MD,MD), F(MD), F2(MD), 1 PHIB(MAXZ,MAXR), PHIP(MAXZ,MAXR), PHIO(MAXZ,MAXR) C DO 2 0 I = 1, M DUM = U(I) DELU = l.D-8 * U(I) IF ( DELU .EQ. 0. ) DELU = l.D-8 U(I) = U(I) + DELU CALL CYLINDER(U,ILAY,IRAD,ISEG_B,ISEG_P,F2,PHIB,PHIP,PHIO) DO 10 J = 1, M DFDU(J,I) = (F2(J)-F(J)) /DELU 10 CONTINUE U(I) = DUM 2 0 CONTINUE C RETURN END C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * DOUBLE PRECISION FUNCTION SSQ(F,M) C Evaluates the sum of squares of the f u n c t i o n s . *************************************************************** IMPLICIT REAL*8(A-H,0-Z) PARAMETER (MD=2000) DIMENSION F(MD) SUM = O.DO DO 10 I = 1, M SUM = SUM + F(I)**2 10 CONTINUE SSQ = SUM RETURN END C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 137 SUBROUTINE INITIALIZE(M,ILAY,IRAD,U) I n i t i a l i z e s the p o t e n t i a l values, sets step p r o b a b i l i t i e s , sets indexing values, and c a l c u l a t e s l a t t i c e s i t e s i n each r a d i a l . * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * IMPLICIT REAL*8(A-H,O-Z) REAL*8 L,L1,L0,LR,LR0,LR1 PARAMETER (MD=2 000,MAXR=2 0) DIMENSION U(MD), L(MAXR), LR(MAXR), LRO(MAXR), LR1(MAXR) COMMON/INDEX/ISTAR1,ISTAR2,ISTAR3 COMMON/SITES/PI,L COMMON/STEPPROB/L1,LO,LR,LRO,LR1 DO 10 I = 1, M U(I) = 0.D0 CONTINUE DO 15 I = 1, IRAD U(IRAD*ILAY+I) = l.D-1 CONTINUE Assign the l a t t i c e step weighting f a c t o r s . LI = 2.5D-1 LO = 5.D-1 DO 20 1 = 1 , IRAD LR(I) = (DBLE(I)-1.DO) / (4.D0*DBLE(I)-2.DO) LR1(I) = DBLE(I) / (4.D0*DBLE(I)-2.DO) LRO(I) = 1.D0 - LR(I) - LR1(I) CONTINUE LRl(IRAD) = LR(IRAD) LRO(IRAD) = 1.D0 - 2.D0*LR(IRAD) C a l c u l a t e the indexing p o s i t i o n s , and the number of s i t e s per r a d i a l i n the c y l i n d e r . ISTAR1 = 0 ISTAR2 = ILAY * IRAD ISTAR3 = 2 * ISTAR2 PI = 4.DO * DATAN(1.DO) DO 30 J = 1, IRAD L(J) = PI * ( 2.DO * DBLE(J) - 1.DO ) CONTINUE RETURN END * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * SUBROUTINE CYLINDER(U,ILAY,IRAD,ISEG_B,ISEG_P, F, PHIB, PHIP,PHIO) Subroutine to c a l c u l a t e the volume f r a c t i o n d i s t r i b u t i o n of a polymer brush, g r a f t e d p a r t i c l e and solvent i n a c y l i n d r i c a l geometry, based on the SF model. U and F are indexed vectors representing a 2D d i s t r i b u t i o n . They are organized by l a y e r , such that the f i r s t IRAD values represent the solvent p o t e n t i a l i n each r a d i a l f o r the f i r s t l a y e r , followed by the second l a y e r 138 C f o r IRAD values and so on. Once the l a t t i c e i s completed f o r C the solvent, the sequence i s repeated f o r the brush and then C the p a r t i c l e . The indexing i s s t a r t e d i n the c o r r e c t p o s i t o n C by the ISTART value. * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * IMPLICIT REAL*8(A-H,O-Z) PARAMETER(MD=2000,MAXZ=50,MAXR=20,MAXS=500) REAL*8 L1,L0,LR,LR0,LR1 DIMENSION U(MD), F(MD), PHIB(MAXZ,MAXR), PHIO(MAXZ,MAXR), 1 PHIP(MAXZ,MAXR), LR(MAXR), LRO(MAXR), LR1(MAXR), 2 GEND(MAXZ,MAXR,MAXS), GTERM(MAXZ,MAXR,MAXS), 3 UOHARD(MAXZ,MAXR), UBHARD(MAXZ,MAXR), UPHARD(MAXZ,MAXR) COMMON/PARAMS/CHIBO,CHIBP,CHIBS,CHIPO,CHIPS,CHIOS,BRUSH, 1 PARTICL,NSEGS COMMON/POSITION/IPOS1,IPOS2,IPOS COMMON/INDEX/ISTAR1,ISTAR2,ISTAR3 COMMON/STEPPROB/L1,LO,LR,LRO,LR1 COMMON/BULK/BULKO,BULKB,BULKP C C C a l c u l a t e the solvent volume f r a c t i o n d i s t r i b u t i o n . C ( p h i ( z , r ) = G(z,r) f o r monomeric solvent ) C BULKO = 1.DO DO 30 1 = 1 , ILAY IMO = 1 - 1 DO 20 J = 1, IRAD PHIO(I,J) = EXP ( -U(IMO*IRAD+J) ) 2 0 CONTINUE 3 0 CONTINUE C C Assign the c o r r e c t free segment weighting f a c t o r s f o r the C brush. C ISTART = ISTAR2 ! Index proper l o c a t i o n i n l i s t . DO 50 I = 1, ILAY IMO = 1 - 1 DO 40 J = 1, IRAD GEND(I,J,1) = EXP( -U(ISTART + IMO*IRAD + J) ) 4 0 CONTINUE 5 0 CONTINUE C C C a l c u l a t e the chain end d i s t r i b u t i o n functions f o r f r e e chains. C CALL ENDDIST(ILAY,IRAD,ISEG_B,GEND) C C Assign the c o r r e c t free segment weighting f a c t o r s f o r C t e r m i n a l l y attached chains. This " g r a f t s " the chains to the C f i r s t l a y e r of the l a t t i c e . C DO 70 I = 2, ILAY DO 60 J = 1, IRAD GTERM(I,J,1) = 0.D0 60 CONTINUE 70 CONTINUE 139 DO 75 J = 1, IRAD GTERM(1,J,1) = GEND(1,J,1) 75 CONTINUE C C C a l c u l a t e the chain end d i s t r i b u t i o n f u n c t i o n s f o r g r a f t e d C chains. Combine the chain end weighting f a c t o r s from free C and g r a f t e d chains to r e t u r n the segment d e n s i t y d i s t r i b u t i o n . C CALL VOLFRAC(ILAY,IRAD,ISEG_B,GEND,GTERM,BRUSH, PHIB,BULKB) C C Assi g n the c o r r e c t free segment weighting f a c t o r s f o r the C p a r t i c l e . C ISTART = ISTAR3 DO 150 1 = 1, ILAY IMO = 1 - 1 DO 140 J = 1, IRAD GEND(I,J,1) = EXP( 14 0 CONTINUE 150 CONTINUE C C C a l c u l a t e the chain end d i s t r i b u t i o n f u n c t i o n f o r the C p a r t i c l e . C CALL ENDDIST(ILAY,IRAD,ISEG_P,GEND) C DO 17 0 I = 1, ILAY DO 160 J = 1, IRAD GTERM(I,J,1) = O.DO 160 CONTINUE 170 CONTINUE C C Assign the c o r r e c t segment weighting f a c t o r s f o r the g r a f t e d C p a r t i c l e . i . e . s p e c i f y the p o s i t i o n of the p a r t i c l e C GTERM(IPOS,1,1) = GEND(IPOS,1,1) C C C a l c u l a t e the chain end d i s t r i b u t i o n f u n c t i o n s f o r a g r a f t e d C p a r t i c l e . Combine the chain end weighting f a c t o r s from f r e e C and g r a f t e d p a r t i c l e to r e t u r n the segment d e n s i t y C d i s t r i b u t i o n . C CALL VOLFRAC(ILAY,IRAD,ISEG_P,GEND,GTERM,PARTICL,PHIP,BULKP) C C C a l c u l a t e the i n t e r a c t i o n energies, and the average C hard core p o t e n t i a l f o r the solvent, brush and p a r t i c l e . C ISTART = ISTAR1 CALL HCPOT(U,ISTART,ILAY,IRAD,PHIB,PHIP,CHIBO,CHIPO, CHIOS, 1 UOHARD) ISTART = ISTAR2 CALL HCPOT(U,ISTART,ILAY,IRAD,PHIO,PHIP, CHIBO, CHIBP, CHIBS, 1 UBHARD) ISTART = ISTAR3 ! Index proper l o c a t i o n i n l i s t . -U(ISTART + IMO*IRAD + J) ) 140 CALL HCPOT(U,ISTART,ILAY,IRAD,PHIO,PHIB,CHIPO,CHIBP, CHIPS, 1 UPHARD) C C C a l c u l a t e the f u n c t i o n values based on the segment d e n s i t y C d i s t r i b u t i o n and average hard core p o t e n t i a l s . C DO 270 1 = 1 , ILAY IMO = 1 - 1 DO 260 J = 1, IRAD AVG = ( UOHARD ( I , J) + UBHARD ( I , J) + UPHARD ( I , J) ) / 1 DBLE(NSEGS) F(IMO*IRAD+J) = 1.D0 - 1.DO/(PHIO(I,J)+PHIB(I,J)+PHIP(I,J)) 1 + AVG - UOHARD(I,J) F(ISTAR2+IMO*IRAD+J) = 1.D0 - 1.DO/(PHIO(I,J)+PHIB(I,J)+ 1 PHIP(I,J)) + AVG - UBHARD(I, J) F(ISTAR3+IMO*IRAD+J) = 1.D0 - 1.DO/(PHIO(I,J)+PHIB(I,J)+ 1 PHIP(I,J)) + AVG - UPHARD(I,J) 260 CONTINUE 270 CONTINUE RETURN END C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * SUBROUTINE ENDDIST(ILAY,IRAD,ISEG,GEND) C Subroutine to c a l c u l a t e the chain end d i s t r i b u t i o n f u n c t i o n . C KEY values handle the r e f l e c t i n g boundary c o n d i t i o n s at IRAD C and ILAY, as w e l l as s o l i d boundary at z = 0. * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * IMPLICIT REAL*8(A-H,0-Z) REAL*8 L I , L0, LR, LRO, LR1, KEY1, KEY2A, KEY2B, KEY3A, KEY3B PARAMETER (MAXZ=50,MAXR=20,MAXS=500) DIMENSION GEND(MAXZ,MAXR,MAXS), LR(MAXR), LRO(MAXR), LR1(MAXR) COMMON/STEPPROB/L1,L0,LR,LRO,LR1 C ISMO = ISEG - 1 IF ( ISMO .EQ. 0 ) RETURN DO 90 N = 1,ISMO NPO = N + 1 KEY1 = O.DO KEY2A = O.DO KEY2B = 1.D0 KEY3A = O.DO KEY3B = 1.D0 DO 70 I = 1, ILAY IPO = 1 + 1 IMO = 1 - 1 IF (I .EQ. ILAY) THEN KEY2A = 1.D0 KEY2B = O.DO END IF DO 60 J = 1, IRAD JPO = J + 1 JMO = J - 1 IF (J .EQ. IRAD) THEN 141 KEY3A = l.DO KEY3B = O.DO END IF GEND(I,J,NPO) = GEND(I,J,1) * ((LI * KEY1 * (LR(J) * 1 GEND(IMO,JMO,N) + (LRO(J) + KEY3A * LR1(J)) * 2 GEND(IMO,J,N) + LR1(J) * KEY3B * GEND(IMO,JPO,N))) + 3 ((LO + KEY2A * LI) * (LR(J) * GEND(I,JMO,N) + (LRO(J) 4 + KEY3A * LR1(J)) * GEND(I,J,N) + LR1(J) * KEY3B * 5 GEND(I,JPO,N))) + (LI * KEY2B * (LR(J) * GEND(IPO,JMO,N) 6 + (LRO(J) + KEY3A * LR1(J)) * GEND(IPO,J,N) + LR1(J) * 7 KEY3B * GEND(IPO,JPO,N)))) CONTINUE KEY1 = l.DO KEY3A = O.DO KEY3B = l.DO CONTINUE CONTINUE RETURN END * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * SUBROUTINE VOLFRAC(ILAY,IRAD,ISEG,GEND,GTERM,RNUM,PHI,BULK) Subroutine t o c a l c u l a t e the chain end d i s t r i b u t i o n f u n c t i o n of g r a f t e d chains. Composition Law i s then used to generate the segment d e n s i t y d i s t r i b u t i o n . **************************************************************** IMPLICIT REAL*8(A-H,0-Z) REAL*8 L, L I , LO, LR, LRO, LR1, KEY1, KEY2A, KEY2B, KEY3A, KEY3B PARAMETER (MAXZ=50,MAXR=20,MAXS=500) DIMENSION GEND(MAXZ,MAXR,MAXS), GTERM(MAXZ,MAXR,MAXS), LR(MAXR), 1 LRO(MAXR), LR1(MAXR), PHI(MAXZ,MAXR), L(MAXR) COMMON/STEPPROB/L1,LO,LR,LRO,LR1 COMMON/SITES/PI,L ISMO = ISEG - 1 IF ( ISMO .NE. 0 ) THEN DO 150 N = l,ISMO NPO = N + 1 KEY1 = O.DO KEY2A = O.DO KEY2B = l.DO KEY3A = O.DO KEY3B = l.DO DO 140 1 = 1 , ILAY IPO = 1 + 1 IMO = 1 - 1 IF (I .EQ. ILAY) THEN KEY2A = l.DO KEY2B = O.DO END IF DO 130 J = 1, IRAD JPO = J + 1 JMO = J - 1 IF (J .EQ. IRAD) THEN 142 KEY3A = l.DO KEY3B = O.DO END IF GTERM(I,J,NPO) = GEND(I,J,1) * ((LI * KEY1 * (LR(J) * 1 GTERM(IMO,JMO,N) + (LRO(J) + KEY3A * LR1(J)) * 2 GTERM(IMO,J,N) + LR1(J) * KEY3B * GTERM(IMO,JPO,N))) + 3 ((LO + KEY2A * LI) * (LR(J) * GTERM(I,JMO,N) + (LRO(J) 4 + KEY3A * LR1(J)) * GTERM(I,J,N) + LR1(J) * KEY3B * 5 GTERM(I,JPO,N))) + (LI * KEY2B * (LR(J) * GTERM(IPO,JMO,N) 6 + (LRO(J) + KEY3A * LR1(J)) * GTERM(IPO,J,N) + LR1(J) * 7 KEY3B * GTERM(IPO,JPO,N)))) 13 0 CONTINUE KEY1 = l.DO KEY3A = O.DO KEY3B = l.DO 14 0 CONTINUE 15 0 CONTINUE END IF C C C a l c u l a t e the n o r m a l i z a t i o n constant f o r g r a f t e d molecules. C SUM = O.DO DO 170 J = 1, IRAD JMO = J - 1 DO 16 0 I = 1, ILAY SUM = SUM + L(J) * GTERM(I,J,ISEG) 16 0 CONTINUE 17 0 CONTINUE GR1 = SUM C = RNUM / GR1 BULK = ISEG * C C C C a l c u l a t e the volume f r a c t i o n p r o f i l e of the brush. C DO 200 1 = 1 , ILAY DO 190 J = 1, IRAD SUM = O.DO DO 180 N = 1, ISEG SUM = SUM + GTERM(I,J,N) * GEND(I,J,ISEG-N+1) / 1 GEND(I,J,1) 18 0 CONTINUE PHI(I,J) = C * SUM 190 CONTINUE 2 00 CONTINUE C RETURN END C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * SUBROUTINE HCPOT(U,ISTART,ILAY,IRAD,PHI1,PHI2,CHI1,CHI2,CHIS, 1 UHARD) C Subroutine to c a l c u l a t e the hard core p o t e n t i a l f o r a given C segment type. The i n t e r a c t i o n p o t e n t i a l f o r i s subtracted C from the t o t a l p o t e n t i a l , u' = u - u ( i n t ) 143 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * IMPLICIT REAL*8(A-H,O-Z) REAL*8 L I , LO, LR, LRO, LR1, KEY1, KEY2A, KEY2B, KEY3A, KEY3B PARAMETER (MD=2000,MAXZ=50,MAXR=20) DIMENSION U(MD), UHARD(MAXZ,MAXR), LR(MAXR), LRO(MAXR), 1 LR1(MAXR), PHI1(MAXZ,MAXR), PHI2(MAXZ,MAXR) COMMON/STEPPROB/L1,LO,LR,LRO,LR1 C KEY1 = O.DO KEY2A = O.DO KEY2B = l.DO KEY3A = O.DO KEY3B = l.DO DO 230 1 = 1 , ILAY IMO = 1 - 1 IF (I .EQ. ILAY) THEN KEY2A = l.DO KEY2B = O.DO END IF DO 220 J = 1, IRAD IF (J .EQ. IRAD) THEN KEY3A = l.DO KEY3B = O.DO END IF PHIAV1 = AVGPHI(I,J,PHI1,KEY1,KEY2A,KEY2B,KEY3A,KEY3B) PHIAV2 = AVGPHI(I,J,PHI2,KEY1,KEY2A,KEY2B,KEY3A,KEY3B) UINT = CHI1 * PHIAV1 + CHI2 * PHIAV2 IF (I .EQ. 1) UINT = UINT + CHIS * LI UHARD(I,J) = U(ISTART+IMO*IRAD+J) - UINT 22 0 CONTINUE KEY1 = l.DO KEY3A = O.DO KEY3B = l.DO 23 0 CONTINUE RETURN END C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * SUBROUTINE ENERGY(U,ILAY,IRAD,ISEG_B,ISEG_P,PHIO,PHIB, 1 PHIP,DELS,DELU,DELA) C This subroutine c a l c u l a t e s the entropy change and Helmholtz C energy change based on the reference s t a t e of unmixed C components. The equations are derived from the s t a t i s t i c a l C mechanics of the l a t t i c e model. ("Polymers at Interfaces") ********************************************************************* IMPLICIT REAL*8(A-H,O-Z) PARAMETER(MD=2000,MAXZ=50,MAXR=20) REAL*8 L, L I , L0, LR, LRO, LR1, KEY1, KEY2A, KEY2B, KEY3A, KEY3B COMMON/PARAMS/CHIBO,CHIBP,CHIBS,CHIPO,CHIPS,CHIOS,BRUSH, 1 PARTICL,NSEGS COMMON/POSITION/IPOS1,IPOS2,IPOS COMMON/INDEX/ISTAR1,ISTAR2,ISTAR3 DIMENSION U(MD), PHIB(MAXZ,MAXR), PHIO(MAXZ,MAXR), L(MAXR), 1 PHIP(MAXZ,MAXR), LR(MAXR), LRO(MAXR), LR1(MAXR) 144 C0MM0N/STEPPR0B/L1,LO,LR,LRO,LR1 COMMON/SITES/PI,L COMMON/BULK/BULKO,BULKB,BULKP C C C a l c u l a t i o n of the entropy part of the p a r t i t i o n f u n c t i o n . C SUM = O.DO ISTART = ISTAR1 DO 30 I = 1, ILAY IMO = 1 - 1 DO 20 J = 1, IRAD SUM = SUM + L(J)*PHIO(I,J)*( LOG(BULKO) - 1 U(ISTART+IMO*IRAD+J) ) 20 CONTINUE 30 CONTINUE IF ( BULKB .NE. O.DO ) THEN ISTART = ISTAR2 DO 50 I = 1, ILAY IMO = 1 - 1 DO 40 J = 1, IRAD SUM = SUM + L( J ) * P H I B ( I , J ) * ( LOG(BULKB)/ISEG_B - 1 U(ISTART+IMO*IRAD+J) ) 4 0 CONTINUE 50 CONTINUE END IF IF ( BULKP .NE. O.DO ) THEN ISTART = ISTAR3 DO 70 I = 1, ILAY IMO = 1 - 1 DO 60 J = 1, IRAD SUM = SUM + L ( J ) * P H I P ( I , J ) * ( LOG(BULKP)/ISEG_P - 1 U(ISTART+IMO*IRAD+J) ) 60 CONTINUE 70 CONTINUE END IF DELS = SUM C C C a l c u l a t i o n of the energy part of the p a r t i t i o n f u n c t i o n . C TERM3 = O.DO KEY1 = O.DO KEY2A = O.DO KEY2B = l.DO KEY3A = O.DO KEY3B = l.DO DO 110 1 = 1 , ILAY IF (I .EQ. ILAY) THEN KEY2A = l.DO KEY2B = O.DO END IF DO 100 J = 1, IRAD IF (J .EQ. IRAD) THEN KEY3A = l.DO KEY3B = O.DO 145 END IF PHIOAV = AVGPHI(I,J,PHIO,KEY1,KEY2A,KEY2B,KEY3A,KEY3B) PHIBAV = AVGPHI(I,J,PHIB,KEY1,KEY2A,KEY2B,KEY3A,KEY3B) PHIPAV = AVGPHI(I,J,PHIP,KEY1,KEY2A,KEY2B,KEY3A,KEY3B) SUM = CHIBO * PHIBAV SUM = SUM + CHIPO * PHIPAV IF ( I .EQ. 1 ) SUM = SUM + 2.DO * CHIOS * LI TERM3 = TERM3 + L(J) * PHIO(I,J) * SUM SUM = CHIBO * PHIOAV SUM = SUM + CHIBP * PHIPAV IF ( I .EQ. 1 ) SUM = SUM + 2.DO * CHIBS * LI TERM3 = TERM3 + L(J) * PHIB(I,J) * SUM SUM = CHIPO * PHIOAV SUM = SUM + CHIBP * PHIBAV IF ( I .EQ. 1 ) SUM = SUM + 2.DO * CHIPS * LI TERM3 = TERM3 + L(J) * PHIP(I,J) * SUM CONTINUE KEY1 = l.DO KEY3A = O.DO KEY3B = l.DO CONTINUE DELU = TERM3 DELA = DELS + DELU RETURN END ************************************************************** DOUBLE PRECISION FUNCTION AVGPHI(I,J,PHI,KEY1,KEY2A,KEY2B, 1 KEY3A,KEY3B) C a l c u l a t e s the l a y e r average volume f r a c t i o n f o r use i n energy c a l c u l a t i o n s . Boundary co n d i t i o n s are handled by the KEY values again. ************************************************* IMPLICIT REAL*8(A-H,0-Z) PARAMETER(MAXZ = 50,MAXR=2 0) REAL*8 L I , LO, LR, LRO, LR1, KEY1, KEY2A, KEY2B, KEY3A, KEY3B DIMENSION PHI(MAXZ,MAXR), LR(MAXR), LRO(MAXR), LR1(MAXR) COMMON/STEPPROB/L1,LO,LR,LRO,LR1 IPO = 1 + 1 IMO = 1 - 1 JPO = J + 1 JMO = J - 1 AVGPHI = ((LI * KEY1 * (LR(J) * 1 PHI(IMO,JMO) + (LRO(J) + KEY3A * LR1(J)) * 2 PHI(IMO,J) + LR1(J) * KEY3B * PHI(IMO,JPO))) + 3 ((LO + KEY2A * LI) * (LR(J) * PHI(I,JMO) + (LRO(J) 4 + KEY3A * L R l ( J ) ) * PHI(I,J) + LR1(J) * KEY3B * 5 PHI(I,JPO))) + (LI * KEY2B * (LR(J) * PHI(IPO,JMO) 146 6 + (LRO(J) + KEY3A * LR1(J)) * PHI(IPO,J) + LR1(J) * 7 KEY3B * PHI(IPO,JPO)))) RETURN END * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * SUBROUTINE OUTPUT(ILAY,IRAD,ISEG_B,ISEG_P,U,PHIB,PHIP,PHIO, 1 I_MIN,I_SEC,ITVAR,ENTR,ENGY,TOTENG) Volume f r a c t i o n s , free segment p o t e n t i a l s , and energy changes based on p a r t i c l e p o s i t i o n are p r i n t e d to an output f i l e . The output format f a c i l i t a t e s 3D p l o t t i n g u s i n g O r i g i n . **************************************************************** IMPLICIT REAL*8(A-H,O-Z) PARAMETER (MD=2000,MAXZ=50,MAXR=20) REAL*8 L DIMENSION U(MD), PHIB(MAXZ,MAXR), PHIP(MAXZ,MAXR), L(MAXR), 1 PHIO(MAXZ,MAXR), ENTR(MAXZ), ENGY(MAXZ), TOTENG(MAXZ) COMMON/PARAMS/CHIBO,CHIBP,CHIBS,CHIPO,CHIPS, CHIOS , BRUSH, 1 PARTICL,NSEGS COMMON/POSITION/IPOS1,IPOS2,IPOS COMMON/INDEX/ISTAR1,ISTAR2,ISTAR3 COMMON/SITES/PI,L COMMON/BULK/BULKO,BULKB,BULKP WRITE(14,*) WRITE(14,'(2X, A39,I4 )') 1 'VALUE FOR GRAFTED CHAIN LENGTH WRITE(14, ' (2X, A3 9,14) ') 1 'VALUE FOR PARTICLE CHAIN LENGTH WRITE(14, ' (2X, A39,14) ') 1 'VALUE FOR CYLINDRICAL LATTICE LENGTH WRITE(14,'(2X,A39,14)') 1 'VALUE FOR CYLINDER RADIUS WRITE(14,'(2X,A39,14)') 1 'START POSITION FOR PARTICLE WRITE(14, ' (2X,A39,14) ' ) 1 'END POSITION FOR PARTICLE WRITE(14, ' (2X, A39,14) ') 1 'CENTER POSITION OF CENTER OF PARTICLE WRITE(14,'(2X,A39,G12.4) 1) 1 'BRUSH-SOLVENT INTERACTION PARAMETER WRITE(14,'(2X,A39,G12.4)') 1 'BRUSH-SURFACE INTERACTION PARAMETER WRITE(14,'(2X,A39,G12.4)') 1 'BRUSH-PARTICLE INTERACTION PARAMETER WRITE(14,'(2X,A39,G12.4)') 1 'PARTICLE-SOLVENT INTERACTION PARAMETER WRITE(14, 1(2X,A39,G12.4)') 1 'PARTICLE-SURFACE INTERACTION PARAMETER WRITE(14,'(2X,A3 9,G12.4)') 1 'NUMBER OF BRUSH MOLECULES IN LATTICE WRITE(14, ' (2X,A39,G12.4) ' ) 1 'NUMBER OF PARTICLE MOLECULES , ISEG_B , ISEG_P , ILAY , IRAD , IPOS1 , IPOS2 , IPOS , CHIBO , CHIBS , CHIBP , CHIPO , CHIPS , BRUSH , PARTICL 147 WRITE(14,'(2X,A39,I4)') 1 'NUMBER OF ITERATION VARIABLES ITVAR 30 40 50 60 WRITE(14,*) WRITE(14,*) 'BRUSH DISTRIBUTION ' WRITE(14,*) 'LAYER RADIAL VOLFRAC ISTART = ISTAR2 DO 50 I = 1, ILAY IMO = 1 - 1 DO 30 J = IRAD, 1, -1 WRITE(14,60) I, (ABS(J-l-IRAD) ) ,PHIB(I,J) 1 U(ISTART+IMO*IRAD+J) CONTINUE DO 40 J = 1, IRAD WRITE(14,60) I,(J+IRAD),PHIB(I,J), 1 U(ISTART+IMO*IRAD+J) CONTINUE CONTINUE FORMAT(3X,12,5X,12,3X,G15.5,3X,G15.5) WRITE(14,*) POTENTIAL' 70 75 80 WRITE(14,*) 'PARTICLE DISTRIBUTION (VOL FRAC)' WRITE(14,*) 'LAYER RADIAL VOLFRAC ISTART = ISTAR3 DO 80 I = 1, ILAY IMO = 1 - 1 DO 70 J = IRAD, 1, -1 WRITE(14,6 0) I,(ABS(J-l-IRAD)),PHIP(I,J), L U(ISTART+IMO*IRAD+J) CONTINUE DO 75 J = 1, IRAD WRITE(14,60) I,(J+IRAD),PHIP(I,J), 1 U(ISTART+IMO*IRAD+J) CONTINUE CONTINUE WRITE(14,*) POTENTIAL' 90 95 100 WRITE(14,*) 'SOLVENT DISTRIBUTION (VOL FRAC)' WRITE(14,*) 'LAYER RADIAL VOLFRAC ISTART = ISTAR1 DO 100 1 = 1 , ILAY IMO = 1 - 1 DO 90 J = IRAD, 1, -1 WRITE(14,60) I,(ABS(J-l-IRAD)),PHIO(I,J), L U(ISTART+IMO*IRAD+J) CONTINUE DO 95 J = 1, IRAD WRITE(14,60) I,(J+IRAD),PHIO(I,J), L U(ISTART+IMO*IRAD+J) CONTINUE CONTINUE WRITE(14,*) POTENTIAL' WRITE(14,*) 148 WRITE(14,*) ' BULK (SOLV) (BRUSH) (PARTICLE) WRITE(14,110) IPOS,BULKO,BULKB,BULKP C WRITE(14,*) WRITE(6,*) WRITE(14,*) 'POSITION (S-S*)/kT (U-U*)/kT (A-A*)/kT' WRITE(6,*) 'POSITION (S-S*)/kT (U-U*)/kT (A-A*)/kT' DO I = IPOS2, IPOS, -1 WRITE(14,110) I, ENTR(I), ENGY(I), TOTENG(I) WRITE(6,110) I, ENTR(I), ENGY(I), TOTENG(I) END DO WRITE(14,*) WRITE(6,*) WRITE(14,*) 'POSITION dS(int)/kT dU(int)/kT dA(int)/kT' WRITE(6,*) 'POSITION dS(int)/kT dU(int)/kT dA(int)/kT' DO I = IPOS2, IPOS, -1 DENTR = ENTR(I) - ENTR(IPOS2) DENGY = ENGY(I) - ENGY(IPOS2) DTOT = TOTENG(I) - TOTENG(IPOS2) WRITE(14,110) I, DENTR, DENGY, DTOT WRITE(6,110) I, DENTR, DENGY, DTOT END DO 110 FORMAT (17,3(G15.6)) C WRITE(14,120) I_MIN, I_SEC WRITE(6,120) I_MIN, I_SEC 120 FORMAT (' CPU time:',16,1H:,12,9H [min :s] ) C RETURN END C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * SUBROUTINE MATINV(A,N,AINV,IERROR) C M a t r i x i n v e r s i o n subroutine, using p a r t i a l p i v o t i n g . ********************************************************************* IMPLICIT REAL*8(A-H,0-Z) PARAMETER (MD=2000) INTEGER PIVOT DIMENSION A(MD,MD), AINV(MD,MD), B(MD), PIVOT(MD), 1 QUOT(MD,MD) C C WRITE(6,*) 'IN MATINV' NMO = N - 1 IERROR = 0 C C Carry out standard Gauss e l i m i n a t i o n with p a r t i a l p i v o t i n g C DO 50 K = 1, NMO KPO = K + 1 BIGGEST = DABS(A(K,K)) PIVOT(K) = K DO 10 I = KPO, N AIK = DABS(A(I,K)) IF (AIK .GT. BIGGEST) THEN 149 BIGGEST = AIK PIVOT(K) = I END IF 10 CONTINUE IF (PIVOT(K) .NE. K) THEN DO 20 J = K, N TEMP = A(PIVOT(K),J) A(PIVOT(K),J) = A(K,J) A(K,J) = TEMP 2 0 CONTINUE END IF IF (A(K,K) .EQ. O.DO) THEN IERROR = 1 ! ZERO PIVOT FOUND RETURN END IF DO 40 I = KPO, N QUOT(I,K) = A(I,K) / A(K,K) A(I,K) = O.DO DO 30 J = KPO, N A(I,J) = A(I,J) - QUOT(I,K) * A(K,J) 30 CONTINUE 4 0 CONTINUE 5 0 CONTINUE IF (A(N,N) .EQ. O.DO) THEN IERROR = 1 ! ZERO PIVOT FOUND RETURN END IF Carry out the e l i m i n a t i o n process as above, and then back s u b s t i t u t i o n to determine the elements of theinverse matrix. DO 90 L = 1, N DO 60 I = 1, N B(I) = O.DO CONTINUE B(L) = l.DO DO 70 K = 1, NMO IF (PIVOT(K) .NE. K) THEN TEMP = B(PIVOT(K)) B(PIVOT(K)) = B(K) B(K) = TEMP END IF DO 65 I = K+l, N B(I) = B(I) - QUOT(I,K) * B(K) 65 CONTINUE 70 CONTINUE AINV(N,L) = B(N) / A(N,N) DO 80 I = NMO, 1, -1 SUM = O.DO DO 75 J = I+l, N SUM = SUM + A(I,J) * AINV(J,L) 75 CONTINUE AINV(I,L) = (B(I) - SUM) / A(I,I) 80 CONTINUE C C C C 60 150 90 CONTINUE RETURN END C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * REAL*8 FUNCTION F_TIME( ) C This f u n c t i o n returns the CPU time i n seconds. * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * IMPLICIT REAL*8 (A-H,O-Z) INTEGER*2 hour, minute, second, hundredth CALL GETTIM( hour, minute, second, hundredth ) F_TIME = ((DBLE( hour ) * 3600.) + (DBLE( minute) * 60.) + l DBLE( second) + (DBLE( hundredth ) / 100.)) RETURN END 151
Cite
Citation Scheme:
Usage Statistics
Country | Views | Downloads |
---|---|---|
China | 11 | 0 |
United States | 2 | 1 |
Japan | 2 | 0 |
City | Views | Downloads |
---|---|---|
Beijing | 11 | 0 |
Washington | 2 | 0 |
Tokyo | 2 | 0 |
{[{ mDataHeader[type] }]} | {[{ month[type] }]} | {[{ tData[type] }]} |
Share
Share to: