@prefix vivo: . @prefix edm: . @prefix ns0: . @prefix dcterms: . @prefix dc: . @prefix skos: . vivo:departmentOrSchool "Applied Science, Faculty of"@en, "Chemical and Biological Engineering, Department of"@en ; edm:dataProvider "DSpace"@en ; ns0:degreeCampus "UBCV"@en ; dcterms:creator "Steels, Bradley Michael"@en ; dcterms:issued "2009-04-28T20:21:55Z"@en, "1997"@en ; vivo:relatedDegree "Master of Applied Science - MASc"@en ; ns0:degreeGrantor "University of British Columbia"@en ; dcterms: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."""@en ; edm:aggregatedCHO "https://circle.library.ubc.ca/rest/handle/2429/7658?expand=metadata"@en ; dcterms:extent "6736785 bytes"@en ; dc:format "application/pdf"@en ; skos:note "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 x1, 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 (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) vl / 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 , 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 ,)!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 ZZ ( 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,= 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 , ' ( ^ ) 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 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-,(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? c^=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 - jz20). 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 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 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 = 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) < Hc = 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 Oxide^ 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 "@en ; edm:hasType "Thesis/Dissertation"@en ; vivo:dateIssued "1998-05"@en ; edm:isShownAt "10.14288/1.0058581"@en ; dcterms:language "eng"@en ; ns0:degreeDiscipline "Chemical and Biological Engineering"@en ; edm:provider "Vancouver : University of British Columbia Library"@en ; dcterms:publisher "University of British Columbia"@en ; dcterms:rights "For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use."@en ; ns0:scholarLevel "Graduate"@en ; dcterms:title "Analysis of interactions between finite-sized particles and terminally attached polymer using numerical self-consistent-field theory"@en ; dcterms:type "Text"@en ; ns0:identifierURI "http://hdl.handle.net/2429/7658"@en .