Polymer



10       PRISM Module

The PRISM module is an implementation of Curro and Schweizer's polymer reference interaction site model (PRISM) theory of polymers in the liquid state. It is used to calculate the structure and thermodynamics of polymer melts and blends. The chains making up these systems may be homopolymers or copolymers. Properties calculated by the PRISM module include pair distribution functions, scattering functions, blend phase diagrams, and approximate equations of state.


Theory

Introduction

The problem of determining the structure and thermodynamic properties of polymer melts and blends is of obvious technological importance. The miscibility behavior of a polymer blend, if known, can help determine its suitability for a particular application since the physical properties of a blend are to a large degree determined by whether it is miscible or not. Miscibility behavior is also relevant to interfacial strength, insofar as the strength of an interface between two polymers is related to the degree to which the polymers tend to mix. In practical terms, if a blend of known polymers can be found with the desired properties for a specific application, it obviates the need to undergo an expensive synthetic search for a new polymer for that application.

Historically, various lattice-based theories have been used to calculate the thermodynamics of polymer solutions and blends, the Flory-Huggins theory (Flory 1953) being the simplest and best known of these. In all of these theories, one or more interaction parameters (e.g., or g, depending on the specific theory) determine the phase behavior of the system, i.e., the conditions under which the system is in one phase as opposed to two or more phases. Lattice theories by their nature are very coarse-grained, and tend to treat interactions at the statistical segment level as opposed to the atomic level. In these theories, the interaction parameters are generally taken as given, or are determined empirically from experimental data (Qian et al. 1991). These parameters, in principle, are determined by the microscopic interactions among the various components of the system. However, no clear prescription exists for determining interaction parameters from, say, a set of realistic atomistic potential energy functions.

In addition to being coarse-grained, the simplest lattice models are mean-field approaches, meaning that they usually neglect to explicitly account for one or more of the following effects: chain connectivity, density fluctuations, concentration fluctuations, and compressibility. The Polymer Reference Interaction Site Model, or PRISM, theory, a very recent approach to the problem of melts and blends, extends beyond these limitations.

PRISM theory, developed over the past five years by Schweizer and Curro (1989) and Curro and Schweizer (1987) is an off-lattice continuum theory of polymer melts and blends. It is based on the small molecule RISM theory developed by Chandler and Andersen (1972) and formulated for "ring polymers" (in the context of a Feynman path integral representation of delocalized electrons) by Chandler and coworkers (1984). Although in practical calculations, PRISM theory is itself somewhat coarse-grained, dealing explicitly with "sites" rather than atoms, it is, nevertheless, more nearly atomistic than lattice models. PRISM theory also allows for concentration and density fluctuation effects, as well as the effect of chain connectivity, all of which are neglected in the Flory-Huggins approach.

PRISM theory is based on statistical mechanical theories of liquids, the basics of which provide enough material to fill a book (Hansen and McDonald 1986). Therefore, this chapter does not present a comprehensive overview and derivation of RISM and PRISM theory. Instead, the goal here is simply to present the final formulation of the theory, along with enough background from liquid state theory to make it understandable. For more detailed treatments and discussion of subtleties that may be omitted from this presentation, please refer to the primary literature (Chandler and Andersen, 1972; Schweizer and Curro, 1989; Lowden and Chandler, 1973, 1974). If you are interested only in the principal assumptions behind the theory--what is required for a PRISM calculation, and what results are obtained from a PRISM calculation--see the Condensed summary of PRISM theory section.

A caveat for users of PRISM theory

PRISM theory is a technique at the forefront of polymer science. This means, among other things, that the theory itself is still undergoing development. While PRISM theory is very powerful, it represents a different type of tool compared to well-established methods such as RIS theory. The full implications of PRISM theory and its range of applicability are simply not known at this time, nor is the extent to which its predictions are quantitatively reliable in general. It has been extremely successful (Honnell et al. 1991) in predicting the structure of polyethylene and alkane melts. In conjunction with density functional theory, it has been able to quantitatively determine the melting points of polyethylene and polytetrafluoroethylene (McCoy et al. 1991a and 1991b). However, for more complicated melts and for blends, there are simply not enough data yet to demonstrate the extent to which the PRISM approach can yield quantitative predictions. This is partly due to the relative newness of PRISM theory, but also due to the general lack of reliable united-atom potentials and the difficulty in the past of routinely calculating the intramolecular (conformational) structure of long polymer chains. See the Condensed summary of PRISM theory section for more on these issues.

Apart from the issue of the quantitative accuracy of PRISM predictions is the issue of obtaining solutions to the PRISM equations in a reasonable amount of time. The PRISM equations are highly nonlinear. Finding a solution can sometimes be an extremely slow process, especially for relatively complicated systems. At high densities or within the two-phase region of a polymer blend, convergence to a solution might never occur. We have strived to make our implementation of this complicated theory as user-friendly as possible. Even so, if you perform many PRISM calculations, you will no doubt find yourself frustrated at times in your attempts to find solutions. As PRISM theory matures and better solution methods are discovered, this situation will no doubt improve.

Condensed summary of PRISM theory

PRISM is a continuum (as opposed to lattice) theory of polymer melts and blends. Each polymer chain in the system is regarded as being made up of sites that interact with sites on other chains. These are the so-called "interaction sites" in interaction site models (ISM's) such as PRISM.

For example, real polyethylene (Figure 34) is made of a sequence of methylene groups each containing a carbon atom and two hydrogen atoms, whereas "PRISM polyethylene" is a chain of united-atom CH2 sites, each site being a single sphere characterized by a hard core diameter, . In principle, one could construct a PRISM representation of polyethylene such that there is a site corresponding to each atom in the CH2 group, but increasing the number of sites in this manner would cause a significant increase in computation time (among other problems). (The time required to perform a PRISM calculation increases rapidly with the number of types of sites used to describe a polymer.) In practice, it appears that the intramolecular arrangement of backbone atoms is much more important in determining the structure of polyethylene melts than is the slight asphericity of the CH2 groups in the backbone. Thus, PRISM calculations based on united-atom representations of alkanes and polyethylene have been found to give structures in very good agreement with experiment (Honnell et al. 1991). Of course, in more complicated polymers, a greater number of sites may be unavoidable.

Figure 34 Polyethylene and "PRISM Polyethylene"

Fundamentally, PRISM theory is a theory of the structure of polymer systems in the liquid state. This structure is represented (in Fourier space) in the form of partial scattering functions, S(k), or (in real space) intermolecular pair distribution functions, g(r). The latter functions describe the spatial correlations between sites on different polymer chains. However, more than just the structure is obtainable from a PRISM calculation. A wide variety of thermodynamic properties of the system can be calculated using the pair distribution functions. These include blend phase diagrams, effective parameters as measured by neutron scattering, and equations of state. The system may be a simple polymer melt, a multi-component blend, a block copolymer melt, or, in principle, any other polymer system in the liquid state.

Every theory has certain quantities (inputs) that you must specify in order to produce results. In the case of PRISM theory, the two most important inputs are

1. the single chain intramolecular partial structure factors

for each type of chain in the system, and

2. the site-site interaction potentials in the form of hard core diameters, ab, and potential well depths, ab, for the Lennard-Jones tails.

Depending on the accuracy required for a given calculation, the single chain scattering functions (see Eq. 509 in Chapter 12, RIS Module) may be calculated by using simple (e.g., Gaussian or freely-jointed) chain models, or more complicated (e.g., RIS Monte Carlo) methods. The site-site interactions are, in effect, united-atom potentials. Insofar as a PRISM representation of a chain is valid (e.g., not excessively coarse-grained), these potentials should in principle be calculable from a fully atomistic force field. In practice, little effort has been devoted to discovering systematic methods of deriving united-atom potentials in this way, and a great deal of study still needs to be done. Until such a scheme is developed, one can only try to make reasonable estimates for the site-site potentials based, for example, on the van der Waals radii and nonbonded potentials of the atoms making up each site.

Modern liquid state theory

The background of PRISM theory can best be started with a summary of the modern theory of atomic liquids (Hansen and McDonald, 1986). In structure, liquids are intermediate between gases, which are virtually structureless, and crystals, which are highly ordered. In a liquid, unlike a crystal, atoms do not vibrate about fixed lattice sites, but are in a more disordered environment and free to diffuse. Yet, because the atoms in a liquid are relatively closely packed, local structural correlations arise from packing constraints (Figure 35).

Figure 35 Local Packing in an Atomic Liquid

The most common way of representing the structure of an atomic liquid is by means of the pair distribution function, g(r,r'). This function plays a central role in theories of liquids. It specifies the relative probability of finding an atom at position r, given that another atom is at position r'. In an isotropic atomic liquid, g(r,r') is a function of the magnitude |r-r'| only; that is, it has no angular dependence or dependence on the absolute location within the liquid, and is generally written g(r) (see Figure 36 (a)). At very short distances, g(r) vanishes due to the hard core repulsions between pairs of atoms. At intermediate distances, g(r) has a roughly oscillatory shape corresponding to the packing layers immediately surrounding an atom in a liquid. At long distances, g(r) approaches unity, meaning that all correlations between the positions of far-separated atoms are lost. For comparative purposes, pair distribution functions for a gas and a crystal are shown in Figure 36 (b) and (c).

Figure 36 Pair Distribution Functions for an Atomic Liquid, a Gas, and a Crystal

Modern liquid state theory was developed mainly over the period from the 1950's through the early 1970's. Although various different approaches have been taken, all modern theories are based on the Ornstein-Zernike (OZ) equation:

Eq. 348     

In this equation, is the number density of particles, h(r) is defined as g(r) - 1, and c(r) is known as the direct correlation function. By itself, this equation is exact but trivial in that c(r) is defined to be the function that makes the OZ equation true. To give the equation more content, a closure relation (see below) must be used in addition.

The OZ equation may, alternatively, be expressed as follows:

Eq. 349     

where the asterisks denote convolution integrals (integrals of the form shown in Eq. 348). The physical picture implied by the OZ equation is as follows. In a dense system, the total correlation between two atoms consists of two contributions. One is a direct correlation, unmediated by other atoms, and is represented by the first term on the right-hand side. The second is an indirect correlation, in which one of the atoms is correlated with a third atom, which in turn is correlated with other atoms, with the chain of correlations finally culminating in the second atom of the original pair (see Figure 37). (The direct correlations here do not result simply from "bare" atom-atom interactions as given by the interatomic potential. The effective interactions responsible for c(r) are more like--but are not equivalent to--potentials of mean force in that they are "renormalized" by the surrounding liquid.) One reason that c(r) is a useful quantity is that it has a much simpler structure and is generally more short ranged than g(r) (or equivalently h(r)). This makes approximate representations for calculational purposes easier to obtain.

Figure 37 Direct and Indirect Correlations in a Simple Liquid

The OZ equation contains two unknown functions: c(r) and h(r). To solve it requires another equation, a closure, relating h(r) and c(r). A closure is an approximation that allows the OZ equation to be solved. Several standard closures exist (Hansen and McDonald 1986); the one used here is known as the mean spherical approximation, or MSA:

Eq. 350     

Eq. 351     

In Eq. 350 and Eq. 351, is the hard core diameter, or distance of closest approach, between two atoms, u(r) is the potential energy function outside the hard core, kB is the Boltzmann constant, and T is the temperature. Eq. 350 is exact for hard core systems; it states simply that the probability of finding two atoms separated by a distance less than their distance of closest approach is zero. (Recall that h(r) is equal to g(r) - 1, so that when g(r) = 0, h(r) = -1). However, Eq. 351 is approximate.

When the system of interest is a liquid of hard spheres, the potential, u(r), is zero for r > . In this case, the MSA closure is equivalent to the Percus-Yevick, or PY, closure. No real liquid is made up of hard spheres, of course, but it has been found that in dense liquids, the structure is predominantly determined by the harsh repulsive forces experienced by atoms at small separations. Consequently, a hard sphere liquid is often a good representation of a real liquid, and deviations from the hard sphere structure may be treated as small perturbations (Hansen and McDonald 1986).

The RISM approach for molecular liquids

The structure of liquids made up of rigid molecules is more complicated than that of atomic liquids, due to the loss of spherical symmetry when two or more atoms combine to form a molecule. This loss of symmetry implies that the pair distribution function between two molecules takes on an angular dependence; i.e., g(r) becomes g(r,;r','), where and ' are the orientational coordinates of the molecules whose centers of mass lie at r and r', respectively. Explicit treatment of this angular dependence can be extremely difficult from a theoretical point of view.

An alternative to explicit treatment of molecular orientations is to regard each molecule in a liquid as consisting of a number of sites, which may be, but do not have to be, atomic centers. The simplification comes from stipulating that the interactions between sites on different molecules are spherically symmetric in form.

The liquid structure is then described, not by means of a complicated, orientation dependent molecular pair distribution function, but by a set of site-site pair distribution functions, g(r), where the subscripts and label the type of site. This general formalism is known by the name Interaction Site Model, or ISM. (The simplifications arising from the ISM formalism come with a price--the loss of explicit orientational information in the distribution functions.)

Although an interaction site may correspond to a single atom, this is not essential. Depending on the nature of the molecule, it is often useful to have fewer sites than atoms, such that two or more atoms are represented by a single site. Figure 38 illustrates possible ISMs for a few molecules. The most important criterion in choosing an ISM representation for a molecule is that the ISM capture enough of the molecule's structure and represent its interactions with other molecules sufficiently well for the application of interest.

An interaction site model is called a Reference Interaction Site Model, or RISM, when the sites are treated as hard spheres, possibly with the hard sphere diameters as adjustable parameters. (Here, the term "reference" indicates that the hard sphere system may form a reference system for a perturbation expansion.) Chandler and Andersen (1972) formulated an approximate integral equation theory for RISM molecules. It is based on a molecular generalization of the OZ equation and on an MSA-like closure relation. This RISM theory has proven very successful in calculating the structures of molecular liquids. In its overall structure, RISM theory is the same for both small molecule and polymeric liquids. The details of this theory are discussed in the next section.

Figure 38 Possible Interaction Site Models for Nitrogen, Carbon Tetrachloride, and Benzene

Polymer RISM Theory

The generalized Ornstein-Zernike equation in the interaction site picture is (Chandler and Andersen 1972):

Eq. 352     

where H(r) is the matrix of intermolecular site-site pair correlation functions for all possible pairs of site types, C(r) is the corresponding matrix of direct correlation functions, and (r) is the matrix of intramolecular site-site distribution functions.

In Eq. 352, the site and molecular densities have been absorbed into the various matrices. To be precise, the elements of H(r) are abhb(r); the elements of (r) are ~iab(r)for sites of type and in the same chain (i.e., chain i), and zero for and in different chains; and the elements of C(r) are simply c(r). Here, a and b refer to the number densities of sites of type and , and ~i refers to the number density of sites of all types contained in chains of type i.

The precise definition of the intramolecular distribution function, ab(r), is:

Eq. 353     

Here, the expressions j {} and k {} indicate that the subscript j is summed over all sites of type , and that k is summed over all sites of type . Ni is the total number of sites in a chain of type i, which itself contains sites of types and (as well as, perhaps, others). The meaning of the term "site type" is discussed below. The angle brackets indicate an ensemble average--that is, a weighted average over all possible conformations of the molecule. The scheme for weighting conformations may be based, for example, on a crude but computationally simple Gaussian or freely-jointed chain model, or on a more sophisticated treatment such as a rotational isomeric state (RIS) model.

Note that (r) is absent in the OZ equation for atomic liquids (Eq. 348). Its presence in the generalized OZ equation is due to the fact that the intramolecular arrangement of atoms has an effect on the intermolecular correlations and vice versa. Figure 39 illustrates this graphically. For a rigid molecule with fixed bond lengths, the ab(r) elements of the (r) matrix are delta functions or sums of delta functions. For a flexible polymer, the intramolecular distribution functions are more complicated.

Figure 39 Coupling of Intramolecular and Intermolecular Correlations in a Polymeric Liquid

Because the elements in the generalized OZ equation (Eq. 352) are matrices, this "equation" is really a set of coupled equations. The number of distinct equations is m(m + 1) / 2, where m is the rank of the matrices. This rank is equal to the number of distinct site types in the system. In this context, two sites on a chain are considered distinct if their direct correlation functions with sites on other chains differ from each other. That is, sites j and k are distinct if cjl(r) and ckl(r) differ, where l is the index of a third site. Two sites are equivalent (of the same type) if these direct correlation functions (for all l) are the same. Thus, in a monodisperse melt of simple ring polymers, all sites are equivalent; m is equal to 1, and the generalized OZ equation reduces to only a single integral equation, which is easily solvable in the RISM approximation.

In a melt of linear simple chains, the situation is in principle more complicated in that sites near the chain ends are chemically distinct from those in the chain interior. Thus, to be rigorously correct, the chain would have to be treated as consisting of N/2 distinct sites, where N is the degree of polymerization. (The factor of 1/2 comes from the symmetry of the chain.) This would render the OZ equation unsolvable for practical purposes. For long chains, however, it is a good approximation to neglect explicit end effects and treat all sites in the chain as equivalent. This is the essential approximation made by Curro and Schweizer (1987) that makes RISM theory tractable when applied to linear polymers. For chains more complicated than simple chains such as polyethylene, it may be necessary to include more than one site type in a PRISM representation of a polymer. For example, a vinyl polymer may be regarded as being made up of three distinct site types: a methylene site, an alpha-carbon site, and a sidegroup site.

Note that the generalized Ornstein-Zernike equation (Eq. 352) is not specific to single-component systems but is completely general. For example, the system described by Eq. 352 may be a one-component melt, a multi-component blend, a polymer solution, or a block copolymer melt, etc. Differences among these systems manifest themselves in the structure of the H(r), C(r), and (r) matrices, but the OZ equation retains the same form.

As was the case for atomic liquids, a closure relation is needed to solve (in an approximate way) the generalized OZ equation. The mean spherical approximation, MSA, appropriately generalized to the molecular case, is used:

Eq. 354     

The symbols in Eq. 354 are analogous to those in Eq. 351, except that here there is a separate interaction potential and set of correlation functions for each distinct site pair. In the PRISM module, the potential used is a hard core with a Lennard-Jones tail, as shown in Figure 40.

Figure 40 Hard Core Potential with Lennard-Jones Tail

Solving the PRISM equations

The Ornstein-Zernike equation, Eq. 352, along with the closure, Eq. 354, form the set of PRISM equations. Chandler and Andersen (1972) devised a scheme for solving equations such as these, which was subsequently refined for more efficient numerical calculations by Lowden and Chandler (1974). The results are presented below.

First, it is convenient to work in Fourier space, in which the OZ equation is simply an algebraic relation:

Eq. 355     

where:

Eq. 356     

with similar definitions for

and

.

Using these definitions, Chandler and Andersen (1972) derived the functional:

Eq. 357     

Here, a is the number density in the system of sites of type ; b is the density of type sites; Tr denotes the matrix trace operation; det denotes the matrix determinant; and 1 is the unit matrix. The above functional has the property that its functional derivative, with respect to the various direct correlation functions, vanishes within the hard core diameter. That is:

Eq. 358     

(See p.72 of Hansen and McDonald's Theory of Simple Liquids, 1986 for a definition and review of functional differentiation.)

To numerically solve Eq. 358, the direct correlation functions, c(r), are approximated within the hard core by the following power series expansion:

Eq. 359     

Four terms in the above expansion have been found to be sufficient for practical calculations (Lowden and Chandler 1973). Outside the hard core, the direct correlation functions are given by the MSA, Eq. 354. With this approximation, Eq. 358 becomes a set of coupled, nonlinear, algebraic equations for the aj coefficients. Given an initial guess for these coefficients, a solution is found iteratively.

Alternative closures

There is growing evidence from experiment (Gehlsen et al. 1992) and simulation (Deutsch and Binder 1992) that some predictions of PRISM theory when using the MSA closure (Eq. 350 and Eq. 351) are qualitatively incorrect. In particular, the MSA predicts, in disagreement with the studies cited above, that the critical temperature, Tc, for an isotopic blend (displaying a UCST phase diagram) scales as the square root of the degree of polymerization, N (i.e., Tc ~ N1/2). The classical Flory-Huggins theory, as well as the above studies, indicate a linear dependence of Tc on N. Yethiraj and Schweizer (1992) have reported a new closure, which they refer to as the Reference Molecular MSA or R-MMSA:

Eq. 360     

Here, C0 is the direct correlation function of a hard core reference system (one that contains only the hard core interactions, but not the potential tails) as determined using the ordinary MSA (Eq. 350 and Eq. 351). U is the matrix of potential tail functions. Unlike previous closures used for molecular and polymeric systems, the R-MMSA directly incorporates the intramolecular structure as manifested in the Omega matrix. Thus, in this formulation of PRISM theory, the single chain structure appears in both the Ornstein-Zernike equation and the closure.

The Lowden-Chandler method described in the Solving the PRISM equations section cannot be used for the R-MMSA. (However, the Lowden-Chandler scheme is used to solve the hard core reference system.) Instead, the PRISM/R-MMSA equations may be solved using a so-called Picard iteration technique, based on fast Fourier transforms. See Hansen and McDonald (1986, p. 129) for further details on this method of solving integral equations. (For the R-MMSA, the iteration is performed in terms of *C*. This differs slightly from the Picard method used for atomic liquids, which performs the iteration in terms of C.)

Two additional closures now available are the Percus-Yevick (PY) and Reference Molecular Percus-Yevick (R-MPY) closure. The PY closure is used for continuous potentials (e.g., full 6-12 Lennard-Jones) rather than the hard core potentials used in the MSA and R-MMSA. The R-MPY closure is a "molecular" version of the PY closure analogous to the R-MMSA (Yethiraj and Schweizer 1992, 1993; Schweizer and Yethiraj 1993).

Numerical details

This section contains a brief discussion about obtaining numerical solutions to the PRISM/MSA equations using the Lowden-Chandler method. Because of their nonlinearity, these equations are mathematically "poorly behaved" in many ways. A Newton-Raphson root-finding algorithm may be used to obtain the final solution to Eq. 358. However, more robust algorithms are often needed to get the aj coefficients sufficiently close to their final values that the Newton-Raphson method does not fail.

Note that finding a solution to Eq. 358 is equivalent to finding an extremum (in this case, a minimum) of IRISM with respect to the aj coefficients. Accelrys research has found empirically that a downhill simplex minimization scheme applied to IRISM is quite effective in moving the coefficients from their initial guess values to the vicinity of their final values. When the simplex method is used in conjunction with a conjugate-gradient scheme, convergence is sometimes faster. Thus, the procedure used in Accelrys's PRISM implementation is as follows. Starting from an initial guess for the aj coefficients, the procedure alternates between several simplex and several conjugate-gradient minimization steps until the derivatives of IRISM with respect to these coefficients are "reasonably small." Then, Newton-Raphson iterations are performed until these derivatives are effectively zero, and Eq. 358 is thus solved.

The above paragraphs describe the Lowden-Chandler procedure for solving the PRISM equations using the MSA. If the R-MMSA or R-MPY closure is chosen instead, the above method is used to solve the hard core reference system only (that is, a system equivalent to the system of interest, except with all the potential tails set to zero). This provides C0 in Eq. 360 above. Then, Eq. 361 is used in a Picard iteration scheme to solve the PRISM equations. In the case of the PY closure, only a Picard iteration scheme is used to solve the PRISM equations.

After a prospective solution is found by one of the methods outlined above, verification is performed to determine if it is indeed a valid solution by confirming that the argument of the logarithm is positive for all k values in Eq. 357. In some cases, a "solution" is found with a negative argument at k = 0. For a blend, this generally means that the system is inside the spinodal (i.e., in the two-phase region) where the PRISM equations, which implicitly assume a homogeneous, one-phase system, are no longer valid. In some cases, it may instead imply that the system is in a gas-liquid two-phase coexistence region, where again the theory no longer applies. These phenomena are mentioned to emphasize that solutions to the PRISM equations cannot necessarily be found for all possible system conditions.


Implementation

The PRISM module allows you to compute a number of distinct properties of dense polymer systems. In this section, the methods used by the PRISM module to calculate these properties are outlined. See the Methodology section for information on how to go about computing specific properties for a given system.

Structure

As previously mentioned, PRISM theory is fundamentally a theory of structure. Quantities such as partial scattering functions and pair correlation functions are calculated from PRISM results in a straightforward way. The expressions used by the PRISM module to compute these functions are:

Eq. 361     

and:

Eq. 362     

Here,

is the matrix of partial scattering functions. These may be combined into total scattering functions if the scattering form factor for each type of site is known. Scattering functions may be compared directly with the results of (static) neutron, X-ray, or light scattering experiments. The

matrix is Fourier transformed to yield the pair distribution functions, g(r). Pair distribution functions may be compared with results of computer simulations. The purpose of performing such comparisons is to probe the validity of the PRISM theory in a particular context, the validity of the potentials used, and/or the validity of the intramolecular distribution functions used in the PRISM calculation. It is possible (currently by trial-and-error) to work backwards from measured or simulated distribution functions to optimize the parameters used in a PRISM calculation.

Equation of state

The equation of state of a one-component system relates the pressure, density, and temperature to one another. The familiar ideal gas relation, P = kBT, is an example of an equation of state. PRISM theory may be used to numerically calculate approximate equations of state for dense polymer systems. In a multicomponent system, the pressure is not only a function of T and , the temperature and density, but also of composition.

Even in the context of PRISM theory, the equations of state computed by the PRISM module are only approximate because approximations beyond those inherent in PRISM theory are necessary if the pressure is to be calculated from only the pair distribution functions and intermolecular potentials. A variant of a virial pressure equation derived by Honnell et al. (1987) is used in the PRISM module to compute the pressure, P:

Eq. 363     

where chain is the total number density of chains in the systems, a and b are site densities, and u(r) is the interaction potential between sites and sites.

In the above equation, three-body contributions to the pressure are omitted. This is done for several reasons. First, the inclusion of three-body terms in the context of PRISM theory requires a "superposition approximation" because PRISM produces only two-body correlation functions. Such an approximation is likely to be unreliable. Second, in calculations that do include three-body terms, it has been found empirically that these terms make only a small contribution to the total pressure (Schweizer and Curro 1988). Finally, without adjustable parameters, the pressure is one of the most difficult quantities to determine quantitatively from either theory or atomistic simulation. The known inaccuracies in virial pressures determined from liquid state theory do not appear to warrant the inclusion of complicated three-body terms.

RPA chi parameter for blends

Scattering functions obtained from small angle neutron scattering (SANS) experiments, together with the random phase approximation (RPA) introduced by de Gennes (1979), may be used to compute effective parameters for polymer blends. It should be noted that, in general, the RPA parameter is not the same as the thermodynamic (in a free energy expression of the Flory-Huggins type). Because of this, and because PRISM theory as such imposes no assumption of incompressibility on a system, the value of the RPA obtained from a PRISM calculation may be a poor indicator of system miscibility. A direct computation of the spinodal curve (as described in the Spinodal curve for blends section) is a direct, and more reliable, indicator of miscibility.

A derivation of the expression for the RPA parameter, in terms of quantities calculated by means of PRISM theory, is presented by Schweizer and Curro (1989). Only the result, slightly rewritten, is presented here as:

Eq. 364     

where R is the site volume ratio, defined as R ¯a/¯b; ¯a and ¯b are the site volumes, and the quantities are the k = 0 values of the direct correlation functions.

Note that the relation in Eq. 364 contains the direct correlation functions rather than the site-site interaction energies. The latter are the quantities that appear in the "bare" 0, which is the Flory-Huggins generalized to the off-lattice case. Thus, RPA represents an interaction parameter that is renormalized by the presence of other sites in the system. Unlike the classical Flory-Huggins 0 parameter, which is simply proportional to 1/T, RPA, in general, has a more complicated temperature dependence. In addition, RPA may depend on blend composition (i.e., volume fraction) as well as the degree of polymerization.

Spinodal curve for blends

The spinodal curve for a blend system represents the boundary of stability of the one phase (miscible) region of the phase diagram with the two phase (immiscible) region (see Figure 41). For further discussion of the spinodal curve, its distinction from the binodal curve, and phase diagrams in general, see Chapter 8, Phase_Diagram Module.

Figure 41 Spinodal Curve of the UCST Type for a Blend

The spinodal represents a second-order phase transition, and as such is characterized by long wavelength fluctuations in the system. The signature of the spinodal is a zero wave vector (k = 0) divergence in the partial scattering functions given by Eq. 361. This corresponds to the condition (Schweizer and Curro 1989):

Eq. 365     

In the PRISM module, the value of this quantity is referred to as the spinodal criterion. If it is positive, the system is in the one-phase, or possibly metastable, region. The closer the spinodal criterion is to zero, the closer the system is to the limit of one-phase

stability. In the two phase spinodal region, PRISM theory breaks down. Thus, in performing a PRISM calculation of the spinodal curve, this region must be approached from the "outside". That is, it is essential that the initial system conditions be in the one phase region, after which the temperature may be varied to approach the spinodal. Further practical considerations are discussed in Methodology, and examples of calculating spinodal curves are given in the online Tutorial Lesson 4. (Refer to the Pilot interface for this tutorial. For more information, see the Tutorial section.)


Command summary

The PRISM module is used to calculate structural and thermodynamic properties of polymer melts and blends based on the Polymer Reference Interaction Site Model of Schweizer and Curro (1989). The most important quantities needed for such calculations are the intramolecular structure factors for each polymer chain and the potentials describing the interactions between sites on different chains in the system. Among the properties that may be computed using the PRISM module are equations of state, effective chi parameters, and phase diagrams (spinodals).

In addition to the core pulldowns in the top menu bar, the PRISM module adds the PR_System, PR_Compute, Pseudo_Atom, Textfile, Graph, and Background_Job pulldowns to the lower menu bar.

PR_System pulldown

Commands in this pulldown pertain to defining the system on which PRISM calculations are to be performed. The PR_System pulldown contains the Components, Pseudo_Atom, Compute_Omega, Omegas, and Interactions commands.

Components

The Components command is used to define the components of the system on which you wish to perform PRISM calculations. The molecules on the screen are used as visualization aids only; none of their structural information is used in the calculation. The sigma (hard core diameter) and epsilon (Lennard-Jones well depth) values are parameters specifying the interactions between sites on different chains. In the Components command, only energy interactions between pairs of equivalent sites are specified; the Interactions command must be used to specify the sigma and epsilon values for non-equivalent sites.

Compute_Omega

The Compute_Omega command computes the partial intramolecular structure factors (omega functions) for a small molecule or oligomer. These functions can then be used as input to PRISM calculations.

Omegas

The Omegas command is used to define the intramolecular structure factor function) for each of the components in the system. This can be for a Gaussian, Freely_Jointed, or Other type of chain (most likely an RIS chain in the last case). If you choose one of the first two, you must also specify the segment lengths for bonds between like sites. An arithmetic combination rule is used for segments between unlike sites:

Eq. 366     

Interactions

The Interactions command allows you to specify Lennard-Jones tail and hard core interactions between sites of unlike type, either on the same molecule or different molecules.

PR_Compute pulldown

Commands in the PR_Compute pulldown are used to define the specific PRISM calculations to be performed and to run the PRISM background job. This pulldown contains the commands: Process, Conditions, Strategy, and PR_Run.

Process

The Process command specifies what it is you wish to calculate and what process--variation of temperature, volume fraction, degree of polymerization, or a combination--is to be followed in performing these calculations. In a multicomponent system, the volume fractions of at most two components are allowed to vary over the course of a calculation.

Conditions

The Conditions command is used to specify whether a PRISM calculation is to be performed at constant pressure or constant density. Using the init, final, and step parameters, you can perform calculations at a number of different pressures or densities. For constant density calculations, you specify whether or not the density is to be corrected for intramolecular overlap in the intramolecular structure factors ( functions). If the Compute Mode selected in the Process command is Spinodal, only a single pressure may be used for constant pressure calculations. In a constant density calculation on a multicomponent system, the system density is computed from the melt densities and volume fractions of the individual components.

Strategy

The Strategy command specifies the strategy to be used to solve the PRISM equations. This includes the choice of closure relation, of minimizers, and of convergence criteria. It is generally only necessary to invoke the Strategy command when you wish to change closures or are having trouble obtaining convergence in a PRISM calculation. Reasonable default values of the various parameters are provided, but no single strategy is optimum for all situations.

PR_Run

The PR_Run command invokes the PRISM background job.

Pseudo_Atom pulldown

The Pseudo_Atom pulldown contains commands to create and modify pseudoatoms. It contains the following commands: Define, Rename, Delete, and List.

Pseudoatoms are defined as the instantaneous average of the coordinates of a set of real atoms. For example, the centroid of a phenyl group could be displayed as a pseudoatom created from the six carbon atoms in the ring. Pseudoatoms are referenced in the same way as atoms: their names can be explicitly specified at the time of their creation or can be generated automatically.

In the DeCipher and Analysis modules, pseudoatoms can be defined and used to study geometric properties. In the DeCipher module, pseudoatoms can be used to conjunction with the commands in the Functions and Geometrics pulldowns to plot and visualize several user-definable properties. In the Discover module, pseudoatoms can be defined and used to define NOE constraints.

Pseudoatoms can also be used in the Color, Display, Label, and List commands in the Molecule pulldown; the Distance, Angle, Dihedral, and XYZ commands in the Measure pulldown; the Center, Move, Overlay, and Superimpose commands in the Transform pulldown; and the commands in the Subset pulldown.

Please see the Insight II User Guide and online documentation for details.

Textfile pulldown

The Textfile pulldown allows you to view the contents of any file containing only text, such as a .log or .outfile file. Textfile contains one command, Get.

Get

The Textfile/Get command presents a list of files that have file extensions matching the search criteria specified by the module in which you are operating. You can select one of the files or type in the name of any other text file. Get automatically pops the textport, then the UNIX more command is used to display the file. When you finish examining the file, a final <Enter> pushes the textport.

Graph pulldown

The Graph pulldown includes the following commands: Boolean, CharSize, Color, Contour, Correlate, Differentiation, Equation, FFT_Real, Get, Info, Interpolation, Integration, Label, Line_Fit, Modify_Display, Move_Axis, Put, Scale_Axis, Smoothing, Threshold, and Tick_Mark.

Please see Chapter 19, Graph Pulldown for more information.

Background_Job pulldown

The Background_Job pulldown allows you to set up background jobs to run concurrently or interactively with Insight II. You are given the choice of whether to send background jobs to a local or remote host. This pulldown is generic and is found in many Insight II modules that run background jobs. The Background_Job pulldown contains the following commands: Setup_Bkgd_Job, Control_Bkgd_Job, Completion_Status, and Kill_Bkgd_Job.

Please see Chapter 15, Background_Job Pulldown for more information.


Methodology

A large number of options are available for performing a calculation using the PRISM module. You can choose among these options based upon what you wish to calculate, the data (such as intramolecular structure factors or united-atom potentials) that you already have available, and the extent to which you require quantitative, as opposed to simply qualitative, results. This section contains general procedures for choosing among the options to perform various types of PRISM calculations, but by no means are all the possibilities exhausted. Because PRISM theory is as new as it is, many issues (such as the best way to derive site-site potentials) remain unresolved.

General considerations

Before performing a PRISM calculation on a particular system, you must first determine whether PRISM is a suitable method for treating the system. In the current formulation of PRISM theory, the dominant forces must be short-ranged, van der Waals type interactions. If your system is dominated by long-ranged electrostatic or dipolar forces, the PRISM module is probably not suitable. Although short-ranged specific interactions, such as hydrogen bonding, are not explicitly treated by the PRISM module, there are ways of choosing the interactions in a blend such that LCST (Lower Critical Solution Temperature) spinodal curves of the type seen in hydrogen-bonded systems are obtained.

If your system can be adequately described by van der Waals forces, then you must next choose a site representation for the molecules in the system. That is, you must determine a reasonable way of describing each molecule as a sequence of connected sites, where each site will probably represent more than one atom. See the Theory section for examples of this process.

An important point to keep in mind is that the computation time required for a PRISM calculation scales roughly as the fourth power of the total number of site types in the system.

Thus, for initial calculations on a system, it is generally advisable to represent the system with a bare minimum of distinct site types. Having obtained preliminary results, if you then think a more sophisticated site representation is required, go ahead and add more sites. For example, if you wish to perform calculations on a vinyl polymer, you might start with a two-site (or even one-site) representation, and only later proceed to a three-site representation, if needed (see Figure 42).

Having chosen a site representation, you must determine the effective interactions between pairs of sites. The PRISM module treats each such interaction as an impenetrable hard core with diameter ab, attached to which is a Lennard-Jones tail, with well depth ab (see Figure 40). Assuming that the sites in your system each represent more than a single atom, these site-site interactions are united-atom potentials.

Figure 42 Two Possible PRISM Representations of PVC

United atom potentials have been published for many small molecules such as methane, but are relatively rare for the groups making up polymers. Ryckaert and Bellemans (1978), as well as others, have published united-atom potentials appropriate for polyethylene, but little work has been done for more complicated chains. Still less has been done to derive and validate methods of determining united-atom potentials from, e.g., a fully atomistic forcefield or quantum mechanics.

This section contains some general suggestions for estimating site-site potentials. One simple and quick approach is to take as the hard core diameter the distance between the center of a site and the outermost point on the van der Waals surface of all the atoms represented by the site. This approach may be useful for determining the relative hard core diameters of different sites, but is likely to overestimate the diameter for any single site. A more systematic approach is to begin with an atomistic forcefield describing the interactions among two sites; determine a spherically averaged (or, alternatively, minimum) energy as a function of separation between the two sites, and determine a best fit of this function to the hard core/Lennard-Jones representation used by the PRISM module. Finally, if experimental data--scattering, equation of state, or phase diagram--are available for systems containing the molecules of interest, it may be best to work backwards to find interaction parameters that yield PRISM results in agreement with the available data. Then, you may proceed to perform PRISM calculations for conditions where experimental data are lacking.

Having chosen a site representation and potential energy parameters for your system, you must next choose a form for the intramolecular structure factor. Again, several options are available. If you wish to perform a relatively quick calculation or would be satisfied with results of a more qualitative, as opposed to quantitative nature, then a Gaussian or freely-jointed chain representation may be adequate. For more quantitative work, it is possible to use the RIS module to perform RIS Monte Carlo calculations of intramolecular structure factors (here called (k), but referred to as S(k) in the RIS module). These structure factors may then be used as inputs for PRISM calculations.

Omega functions for small molecules and oligomers

If you are interested in performing PRISM calculations on systems containing small molecules or oligomers (e.g., polymer solutions), you can use the Compute_Omega command to calculate the omega functions for such molecules. For a flexible molecule, you may use a molecular dynamics trajectory to produce results averaged over the molecule's accessible conformations.

Before using the Compute_Omega command, you must use the commands in the Pseudo_Atom pulldown to specify the scattering centers (sites) in the molecule. For the purpose of calculating the omega functions, pseudoatoms with identical alphabetic parts of their names are treated as equivalent. For example, pseudoatoms named A, A1, and A23 would be treated as equivalent sites and distinct from pseudoatoms named B2, AB4, or C.

The omega functions are computed with the following formula:

Eq. 367     

where N is the total number of sites defined in the molecule; k is the wavevector; the index i is summed over all sites of type a; the index j is summed over all sites of type b, and rij is the distance between sites i and j. If the molecule is flexible, then the above expression must be averaged over several frames in a trajectory, using the Use_Trajectory option.

Once the omega functions are calculated, the name of the file containing these data is written to the omega list file (.omlist file). To make use of these data in a PRISM calculation, you should then specify this file when you use the Omegas command. The DP that you specify for the small molecule in the Omegas command should be the same as that listed in the .omlist file.

Using the PRISM module to specify the system

Having decided on a PRISM representation of your system based on the considerations described above, your next step is to specify the system parameters in the PRISM module.

The chain or chains in the system are built using Polymerizer. You need to build only one chain of each distinct type. PRISM calculations are done for bulk systems, and the chains you build are used only as visualization aids to help you specify the site-site interaction potentials.

The PR_System/Components command in PRISM is used to specify the number of site types in the chain and the hard core and Lennard-Jones parameters for interactions between equivalent sites, for each polymer in the system. You should use the PR_System/Interactions command to specify interactions between non-equivalent sites.

For a given polymer component, each site type is associated with a distinct atom in the chain. Which atom you use to represent a site is unimportant, so pick the one that you find most convenient. Again, no details of the chemical composition of the molecules on the screen are used in the PRISM calculations. The molecules are used as visualization aids, but you specify explicitly all of the information needed for the PRISM calculation.

The PR_System/Omegas command is used to specify the intramolecular structure factor, for each polymer in the system. The name of the Omega List File (see Appendix C, File Formats) must be specified for chains that are not Gaussian or freely-jointed chains (e.g., RIS chains). The Omega List File lists the names of files containing structure factor data for various temperatures and degrees of polymerization.

The PR_System/Interactions command is used to specify the hard core and Lennard-Jones interaction parameters, for each pair of non-equivalent sites in the system. Rather than specifying each interaction separately, you can tell PRISM to use combining rules [AB = (1/2) (AA + BB) and AB = (AABB)1/2] to determine these parameters.

Calculation options

After defining your system, the commands in the PR_Compute pulldown are used to tell the PRISM module what to calculate and under what system conditions to perform the calculations.

The PR_Compute/Process command is used to select among the three general types of calculations:

For any of these options, the initial and final values and step size for the temperature must be specified. In the spinodal case, the initial temperature must place the system in a one-phase region. The PRISM background job locates the spinodal by finding the temperature at which the spinodal criterion (see the Implementation section) vanishes. For all three options, you can choose to perform calculations over a range of compositions (volume fractions) or degrees of polymerization.

The PR_Compute/Conditions command is used to specify whether calculations are to be performed at "constant density" (for blends, this is really constant site volume) or constant pressure.

The constant pressure calculations are performed by varying the density until the desired (virial) pressure is reached. Constant density calculations are performed directly and, thus, are more rapid. Therefore, unless there are special considerations (e.g., you wish to account for thermal expansion), it is generally best to perform PRISM calculations at constant density.

If you have trouble getting the PRISM calculations to converge in a reasonable amount of time, the PR_Compute/Strategy command can be used to optimize the parameters that determine how the background job goes about solving the RISM equations.

Next, the PR_Compute/PR_Run command is used to execute the PRISM background job.

While the job is running, the Textfile/Get command can be used to monitor its progress by reading the log file (extension .prlog). When the background job successfully finishes, you are left with one or more .prtbl files that may be plotted using the Graph/Get command. See Appendix C, File Formats, for descriptions of the contents of these files.

More on site-site potentials

Choosing suitable site-site potentials is perhaps the most challenging aspect of performing a PRISM calculation. As previously mentioned, much more research needs to be done to determine systematic ways of deriving such potentials from first principles. Below are further guidelines to help you get started in estimating potentials for particular types of systems.

Structure determination

If your main interest is in determining the structure of a one-component nonpolar system, it is generally reasonable to use exclusively hard sphere potentials, with the Lennard-Jones well depths set to zero. First, this reduces the number of parameters in the problem. Second, it is known from studies on small molecule liquids, as well as liquid polyethylene, that by far the dominant intermolecular forces in determining the liquid structure are the hard core repulsions (Hansen and McDonald 1986 and Honnell et al. 1991). Thus, any attractive tails tend to have a relatively minor effect on the structure at a given density. However, such quantities as the pressure are significantly affected by the presence of attractive tails. As mentioned above, you can estimate the hard sphere diameters from the atoms' van der Waals radii or can adjust these to yield agreement with experimental or simulation data (e.g., by matching the location of the first peak in g(r) ).

Spinodal curve determination

If you wish to compute the spinodal curve of a blend, you must include Lennard-Jones tail interactions in the PRISM representation of your system. Otherwise, the system is an athermal blend, whose miscibility depends only on density and volume fraction, but not on temperature. The simplest, although somewhat artificial, way to obtain a UCST phase diagram for a binary (A-B) blend of simple chains is to set AA = BB = 0 and AB < 0. The negative value of AB implies repulsion between A sites and B sites, which tends to favor phase separation at low temperatures. (Recall that is the well depth. Therefore, a negative implies the presence of a repulsive "hill" rather than an attractive well.) Note that in a system such as this, the temperature affects the structure only through the AB term. Thus, if AB is doubled, for example, then points on the phase diagram are simply shifted upward by a factor of two, but no change in the shape of the phase diagram occurs, assuming that the intramolecular structure varies little with temperature This is true for Gaussian or freely-jointed chains, but not for RIS chains.

At this writing, the existence of a spinodal curve of the LCST type arising from differences in thermal expansion coefficients (or free volume differences) of the two system components has not been investigated in the context of PRISM theory. However, there is a way to mimic an LCST arising from "specific interactions" (Coleman et al. 1991) in a binary blend by judicious choice of and parameters. For example, one may set AA = BB < AB, with AA = BB = 0 and AB > 0. If the difference AB - AA is large enough, then the athermal system (corresponding to the high temperature limit, where the potential tails become irrelevant) will phase separate. Similarly, if AB is sufficiently large, then the A-B attraction (mimicking in very simple fashion a specific interaction, such as hydrogen bonding) will cause mixing at low temperature, resulting in LCST behavior. Apart from studies performed at Biosym (now Molecular Simulations) confirming the existence of this phenomenon (Honeycutt 1991), little research has been done in the context of PRISM theory for systems of this type.

If you have trouble finding a solution...

When performing calculations on a new system, the PRISM equations may be slow to converge or may not converge to a solution at all. With the Verbose_Output parameter of the PR_Compute/PR_Run command set to On, you can monitor the progress of a PRISM calculation by repeatedly inspecting the .prlog file using the Textfile/Get command. If convergence is simply slow, you can:

1. use the first solution you obtain as the initial guess for subsequent calculations (by the use of the Use_Previous_Soln parameter of the PR_Compute/PR_Run command) and/or

2. adjust the default optimization parameters using the PR_Compute/Strategy command.

As an example regarding item two, the simplex minimizer and conjugate gradient minimizer often complement each other well. However, each iteration of the latter minimizer is considerably slower than the former, and sometimes the conjugate gradient stage is not needed. To test this, you can set the Conjugate parameter of the Strategy command to off and observe whether subsequent calculations proceed more rapidly.

If you are unable to obtain a solution at all for a blend system, or a one-component system with attractive interactions, then it is possible that the system is in a two-phase region. If you expect your system to exhibit UCST behavior, then try raising the temperature. If you expect LCST behavior, try reducing the temperature. You may of course vary other quantities, such as the density (particularly for athermal blends) or volume fraction, in an attempt to find a solution. In the one component case (and sometimes in the two component case), two phase behavior corresponds to liquid-vapor phase separation. Try raising the temperature or reducing the depth of the attractive wells.

Another possibility, if no solution is found for a one or multicomponent system, is that the density or the hard core diameters are too large. (In some cases, the PRISM background job can determine that this is so and write a corresponding message to the .prlog file.) Try reducing one or more of these quantities.

Further practical considerations

You are reminded again that the time required to perform a PRISM calculation increases drastically with the number of site types in the system. Table 11 shows the dependence of computation time on the number of site types for a Gaussian polymer melt. (For n site types, the system was modelled as an n-component system with the different components having identical parameters. Computations were performed on an SGI 4D/25.) When in doubt, use fewer rather than more sites for exploratory calculations.

Table 11      Computation Time vs. Number of Site Types

n t(min)
1   .3  
2   3.9  
3   19.9  
4   169.0  

For a Gaussian or freely-jointed chain containing more than one type of site, the time required for the PRISM program to calculate the intramolecular structure factors scales as the square of the degree of polymerization. For large chains, this step may take up a large fraction of the total computation time. Because of this, an omega list file is automatically written to disk when structure factors of this type are calculated. The name of the omega list file is <run_name>.omlist. In subsequent calculations, you can then use the Other option for the Omega Type parameter of the PR_System/Omegas command, and specify this file as the Omega List File. For multiple PRISM calculations on large chains, this can save a great deal of time.

When calculating spinodal curves for blend systems with deep Lennard-Jones wells (i.e., large, positive values), it is possible to encounter a liquid-vapor spinodal rather than the desired liquid-liquid spinodal. The PRISM program currently is unable to distinguish between these situations. One signature of a liquid-vapor phase separation is a very flat, or even concave upward, apparent "UCST" spinodal curve (when plotted as a function of volume fraction). By performing a spinodal calculation for a single component system containing only one of the components in the original system (e.g., the component with the deepest Lennard-Jones wells), you can generally confirm whether this is the case, because the only (zero wave vector) spinodal possible in a one-component PRISM system is a liquid-vapor spinodal. Thus, if you get phase separation in the one-component system at a temperature near that at which phase separation is predicted in the blend, the blend spinodal you have obtained is probably a liquid-vapor spinodal. You can more directly confirm the presence of this phenomenon by looking at the compressibility in the blend's .prlog file. At a liquid-vapor transition, the compressibility diverges.

Final comments on closures

You may find yourself having trouble deciding whether to use the MSA or R-MMSA closure in a PRISM calculation. While closure relations for use in PRISM calculations remain the subject of active research, some general guidelines are as follows:

For polymer blends with temperature-dependent interactions (i.e., potential tails), the R-MMSA would appear to be the better closure for obtaining quantitative values of the critical demixing temperature, Tc, from realistic potential energy functions, and thus is more rigorously valid. However, the Picard iteration scheme is less well-behaved than the Lowden-Chandler method used to solve the MSA version of the theory. Also, in spite of its known inadequacies, the MSA is able to reproduce certain experimentally observed trends, such as effects of microstructure on blend miscibility (Honeycutt 1992b). For melts and athermal blends, the MSA appears to be reasonably adequate.

With regard to the relative validity of the two closures now used in the PRISM module, the case of block copolymers is less clear cut than the case of blends. Calculations by Schweizer (1991) indicate that the R-MMSA yields a second-order microphase separation transition (MST), marked by a divergence in the partial scattering functions at nonzero wavevector. Leibler's (1980) seminal mean field treatment of diblocks predicts the same. However, Leibler and others argue that in real diblock systems, fluctuations should destroy the second order transition (Leibler 1980; Fredrickson and Helfand 1987; Brazhovskii 1975), turning it into a first-order transition. This prediction is consistent with the MSA results, which fail to predict a second-order transition (Schweizer, 1991; Honeycutt 1992a). Thus, the MSA may be more appropriate in some ways for block copolymer systems than is the R-MMSA.


Tutorial


Pilot online tutorials

Many tutorials are available online for use with the Pilot interface.

To access the online tutorials for PRISM, click the biplane or mortarboard icon in the Insight interface, or select Online_Tutorials from the Help pulldown.

Then, from the Open Tutorial window, select Polymer Modeling and Property Prediction tutorials, and then select PRISM Module Tutorials from the list of modules. Choose from the list of available lessons:

Lesson 1: Structure of a Polyethylene Melt

Lesson 2: Equation of State Calculation for Freely-Jointed Chain of Polyethylene

Lesson 3: Scattering Functions of a Block Copolymer

Lesson 4: Spinodal Curve for a Simple Blend

Lesson 5: Generating an LCST Spinodal Curve

Lesson 6: Estimating PRISM parameters: One Approach

(Lessons 7, 8, and 9 are not included as online tutorials. They are given below.)

You can access the Open Tutorial window at any time by clicking the Open File button in the lower left corner of the Pilot window.

For a more complete description of Pilot and its use, click the on-screen help button in the Pilot interface or refer to the Introduction in the Insight II manual.

Setting Up for PRISM Lessons 7, 8, and 9

It is suggested that you perform these examples in an empty directory in which you have write access.

Begin by creating a directory called prism_tut and changing to it. To do this enter the following at the prompt in a UNIX window:


>	mkdir prism_tut

>	cd prism_tut
You will need some auxiliary files to perform this tutorial. These files specify single chain structure factors for polyethylene computed using the RIS Monte Carlo technique. Enter the following commands:


>	polymer_tutorials

>	cp $PRISM_EXAMPLES/*.omlist $cwd

>	ln -s $PRISM_EXAMPLES/skgen skgen


Overview of tutorial lessons

The lessons in this tutorial illustrate several of the capabilities of the PRISM module. The topics covered are:

Some of the lessons are fairly time-consuming, and you will probably not want to do them all in one sitting. But, it is suggested that you do them in order, as later lessons build on what is learned in earlier ones.

Lesson 7: Using RIS Monte Carlo structure factors in a PRISM calculation

For PRISM calculations that presume to be quantitative, good single chain structure factors are essential. In previous lessons, you have several times used RIS generated structure factors (referenced by the file PEomega.omlist) in PRISM calculations. In this lesson, you learn how to generate your own RIS omega files for PRISM calculations.

This lesson assumes a basic familiarity with the RIS module. If you are not familiar with this module, see the RIS Tutorials section in Chapter 12. Note: the function referred to as(k) or omega in the PRISM module is known as S(k) or S_of_k in the RIS module.

This lesson uses the following commands, accessed from the specified module.

Polymerizer module

Homopolymer

Polymerize

RIS module

Sequencer

RIS_Monte_Carlo

Graph

S_Run

Scatterers

Get

Range

Put

MC_Run

Viewer module

Session

Unix

Quit

1.   Create a Poly(vinyl chloride) Molecule

If you are continuing from the previous lesson, make sure that all the objects on the screen have been deleted before continuing. Otherwise, start Insight II by typing insightII at the UNIX prompt.

Click the Accelrys icon and select the Polymerizer module and build a 10-mer of poly(vinyl chloride) (PVC) using the Homopolymer/Polymerize command.

2.   Sequence the Molecule

Click the Accelrys icon and select the RIS module. Select the Sequencer/S_Run command, accept all default values for the parameters and select Execute. When the Sequencer job finishes, proceed to the next step.

3.   Define the Scattering Centers

Select the RIS_Monte_Carlo/Scatterers command. Set the focus to the Scatterer Name parameter, and enter CH2. For the Backbone Atom, use the mouse to pick any unsubstituted methylene carbon atom in the chain backbone. Select Execute.

Again set the focus to Scatterer Name, and, this time enter CCl. Set the focus to Backbone Atom and use the mouse to pick any chlorine-substituted carbon atom in the chain backbone. Pick this same atom again for Scatter Atom 1, and the attached chlorine atom for Scatter Atom 2. Select Execute.

You have just defined scattering centers appropriate for a two site type PRISM representation of PVC. The first site type represents a methylene group. The second, a substituted methylene and chlorine group together, with the site center located halfway between the carbon atom and the chlorine atom.

Depending on the application, two site types may or may not be enough to adequately represent this vinyl chain. It is possible that three site types--methylene, substituted methylene, and chlorine--will be required. However, due to the slowness of PRISM calculations for systems containing large numbers of site types, it is generally best to start with fewer sites, adding more only as necessary.

4.   Specify the Small Wavevector Range

Select the RIS_Monte_Carlo/Range command. Range To Set: S_of_k

S_of_k Mode: Order_N_Squared

S_of_k Mode specifies the algorithm used to compute S(k) (known as(k) in the PRISM context).

Here, N is the degree of polymerization. The order N method for calculating S(k) is more approximate than the order N2 method, but is much faster for large chains (containing more than 100 bonds or so). Since the PVC chain used in this lesson is small, the order N2 method is adequate.

Min Argument: 0

Max Argument: 1.5

Number Of Bins: 150

Select Execute.

You will perform the S(k) calculation in two steps, one for small wavevector and one for larger wavevector. The PRISM module requires structure factors out to k 20 Å-1, but because S(k) varies more rapidly at small k, a finer grid is required for smaller k values than for the larger values. In this step, you have specified the small k fine grid for the first calculation.

5.   Run the First RIS Monte Carlo Calculation

Select the RIS_Monte_Carlo/MC_Run command. Enter pvc_smallk for Run Name and the PVC molecule name for Assembly/Molecule. Toggle s_of_k to on and select Execute.

This calculation takes about 1.5 minutes on a Silicon Graphics 4D/25. When it is finished, proceed to the next step.

6.   Specify the Large Wavevector Range

Select the RIS_Monte_Carlo/Range command. Range To Set: S_of_k

S_of_k Mode: Order_N_Squared

Min Argument: 1.5

Max Argument: 20.0

Number Of Bins: 150

Select Execute.

7.   Run the Second RIS Monte Carlo calculation

Select the RIS_Monte_Carlo/MC_Run command. Enter pvc_largek for Run Name and the PVC molecule name for Assembly/Molecule. Toggle S_of_k to on and select Execute.

This calculation takes about 1.5 minutes on a Silicon Graphics 4D/25. When it is finished, proceed to the next step.

8.   Plot and Write Out the Data

Use the Graph/Get command to plot the data from the file pvc_smallk_Sk.tbl. Make three plots on the same graph in the following order:

S_CH2-CH2(k) versus k

S_CCl-CCl(k) versus k

S_CH2-CCl(k) versus k

Then use Graph/Put to write this graph to a file named pvc_T300a_om.tbl. Repeat the procedure in the previous paragraph for the data in the file pvc_largek_Sk.tbl, putting the three plots on a new graph, and using Graph/Put to write the data to a file named pvc_T300b_om.tbl.

9.   Create the Omega List File

Use the Session/Unix command to open a UNIX shell. At the UNIX prompt, enter head pvc_smallk_Sk.tbl to see the first few lines of this file. The sum of the numbers appearing on the lines beginning with the word Scatterer: is the degree of polymerization for PRISM purposes (19 in this case).

In the PRISM module, the DP of a chain is the total number of sites of all types in that chain. Since end groups are not included in this particular scattering calculation, the DP is 19 rather than 21.

Using any text editor (e.g., vi) create a file named pvc.omlist that contains only the following line:

>	19 300.0 pvc_T300a_om.tbl pvc_T300b_om.tbl

Having done this, you may now use the file pvc.omlist as the Omega List File in the Omegas command of the PRISM module. Note that if you wish to perform calculations on PVC for a DP other than 19 or at a temperature other than 300 K, you will have to compute more RIS Monte Carlo scattering functions. However, the PRISM program is capable of interpolating omega data if none are available for the exact DP and temperature of interest. Temperature interpolation tends to be more reliable than DP interpolation. If interpolation has to be done, a warning is included in the PRISM log file. For more details on the formats of omega files and omega list files, see Appendix C, File Formats.

Continue with the next lesson or use the Session/Quit command to exit Insight II.

Lesson 8: Troubleshooting in PRISM calculations

Sometimes a PRISM calculation is very slow to converge to a solution. There are several possible reasons for this: a large number of site types in the system, an initial temperature far from the spinodal temperature (in the case of a spinodal calculation), the system being too dense or within the spinodal such that no PRISM solution exists for the particular system conditions, a poor optimization strategy, and so forth. The examples in this lesson are designed to illustrate some of these situations, and to show you strategies for dealing with them.

This lesson uses the following commands, accessed from the specified module:

PRISM

PR_System

PR_Compute

Textfile

Background_Job

Components

Process

Get

Kill_Bkgd_Job

Omegas

Conditions

PR_Run

Viewer module

Session

Quit

1.   Perform a Calculation on a Too-dense System

If you are continuing from the previous lesson, delete all the objects on the screen by typing delete *. Otherwise, start Insight II by typing insightII at the UNIX prompt.

Build a polyethylene chain as in step 2 of lesson 1.

Using the PR_System/Components command, define a one component, one site type polyethylene system with Sigma11 equal to 4.5 and Epsilon11 equal to 0.

Using the PR_System/Omegas command, specify a Gaussian distribution with a DP Of Chain of 1000. Select Execute.

Select the PR_Compute/Process command. Specify a Structure_and_Chi calculation at 450 K. Select Execute.

Select the PR_Compute/Conditions command. Set the Condition Mode to Const_Density of 0.8. Note that constant density is set to 0.6. Select Execute.

Finally, start the calculation with PR_Compute/PR_Run.

After a few seconds, the following message:


Runaway solution: system is probably too dense or hard cores too large

appears at the bottom of the screen. The problem is that the hard core diameter is too large for a methylene group, resulting in an artificially large packing fraction at a density of 0.8 g cm-3. The solution to the problem is to reduce Sigma11 to a more reasonable number such as 3.9 Å and to repeat the calculation.

2.   Perform a Spinodal Calculation with a Poor Initial Temperature

Set up a blend system as in steps 1 through 4 of Lesson 4 given in Pilot. In the PR_Compute/Process command, set these parameters: Compute Mode: Spinodal

Tinit: 500

Tfinal: 100

Tstep: -20

Phi1init: 0.2

Phi1final: 0.2

Conditions: Const_Density

In the Strategy command: Conjugate: off

Closure: MSA (default value)

Use the PR_Compute/PR_Run command to start the calculation.

Monitor the progress of the calculation by using the Textfile/Get command to read the .prlog file for your job. It takes several minutes to find a PRISM solution at the initial temperature, but calculations proceed more rapidly after that. Note, though, that the SPINODAL criterion (which vanishes at the spinodal) is large, and decreases only slowly as the temperature is dropped. This is because the initial temperature far exceeds the spinodal temperature.

Note also that the SPINODAL criterion does not necessarily vary monotonically with temperature far away from the spinodal. Because of this, the PRISM job may have trouble locating the spinodal if the initial temperature differs too greatly from the spinodal temperature. When you see this occurring, it is often best to kill the background job and restart the calculation with a more suitable initial temperature. Use the Append_Output option in the PR_Compute/PR_Run command, if necessary, to avoid overwriting any spinodal data generated by the initial calculation.

To kill the PRISM background job, select the Background_Job/Kill_Bkgd_Job command. Pick the job from the Running_Jobs value-aid and select Execute.

3.   Attempt a Calculation within a Gas-Liquid Spinodal

When performing PRISM calculations on systems with deep attractive tails, you may encounter situations where no convergence occurs in the PRISM calculation. Sometimes, this indicates liquid-vapor phase separation. You have already encountered this situation in Lesson 6 given online through Pilot for the 60 K well case at high temperatures. When this is the case, the remedy is to perform the calculations at higher density or pressure, or alternatively, to decrease the depth of the potential wells if these are not already known to be accurate.

Set up a polyethylene system as in Lesson 1, but with Epsilon11 equal to 72 K rather than 0. Use the RIS Monte Carlo omegas (Omega List File is PEomega.omlist) in the PR_System/Omegas command.

In the PR_Compute/Process command: Tinit: 450

Tfinal: 250

Tstep: -50

Use the PR_Compute/PR_Run command to start the background job.

Use the Textfile/Get command to monitor its progress.

Note that convergence to a solution occurs very rapidly at the first few temperatures, but that no convergence occurs once the temperature falls to 300 K. This is because the system is within the liquid-vapor two-phase region. In a constant density calculation, liquid-vapor coexistence is encountered when the temperature drops, contrary to the constant pressure case. Liquid-vapor coexistence is only possible for a system with attractive interactions. You can verify this by setting Epsilon11 to 0 and noting that you easily find a solution at 300 K and below for this hard sphere system.

Use the Session/Quit command to exit Insight II.

Lesson 9: Computing the DP dependence of the critical temperature with the R-MMSA closure

In this lesson, you compute the critical temperature as a function of the degree of polymerization for a simple model of a symmetric isotopic blend, using the R-MMSA closure.

Note: This lesson takes about 40 minutes on an SGI Indigo (R3000, 33MHz), and longer on slower machines such as an SGI IRIS 4D/25. As in the previous lesson, less detail is provided here than in other tutorials because it is assumed that you are now familiar with the procedures for setting up PRISM calculations.

This lesson uses the following commands, accessed from the specified module:

PRISM

PR_System

PR_Compute

Graph

Components

Process

Get

Omegas

Conditions

Interactions

Strategy

PR_Run

Viewer module

Session

Quit

1.   Start Insight II

Start Insight II by typing insightII at the UNIX operating system prompt.

2.   Build Two Polyethylene Oligomers

Click the Accelrys icon and select the Polymerizer module. Build two polyethylene oligomers, each with a DP of 6.

3.   Set Up the PRISM Calculation

Click the Accelrys icon and select the PRISM module. Select the PR_System/Components command. Add the two molecules to the system, each with Sigma11 equal to 3.7 Angstroms and Epsilon11 equal to 0.

Select the PR_System/Omegas command. Define the Omega functions for each molecule by accepting all default values.

Select the PR_System/Interactions command. Interaction Param: Epsilon

Use_Combining_Rule: off

Pick the site atoms on the two different molecules, and specify an EpsilonAB value of -0.1. Select Execute. While still in the PR_System/Interactions command: Interaction Param: Sigma

Use_Combining_Rule: on

Select Execute.

Select the PR_Compute/Process command. Compute Mode: Spinodal
Tinit: 300

Tfinal: 100

Tstep: -20

Component 1: POLYETHE

Component 2: POLYETHE0

Phi1init: 0.5

Phi1step: 0

Phi1final: 0.5

DP Mode: Vary_All_DP

DPmin: 100

DPmax: 400

DPstep: 100

Select Execute.

Note that all computations will be performed at a volume fraction of 0.5. For a symmetrical system such as this one, 0.5 is the critical volume fraction. Thus, the spinodal temperature computed by the PRISM module is also the critical temperature.

Select the PR_Compute/Conditions command to specify a Const_Density calculation.

Select the PR_Compute/Strategy command. Toggle Closure to R_MMSA. Leave all other parameters as they are, and select Execute.

4.   Run the PRISM Calculation

Select the PR_Compute/PR_Run command. Enter mmsa_dp_dep as the Run Name, and select Execute.

Note: Unless you have a very fast workstation, this calculation will take at least 30 minutes.

5.   Plot the Critical Temperature vs. DP

When the calculation is done, select the Graph/Get command, and enter mmsa_dp_dep_sp.prtbl as the File Name. X Function: DP

Y Function: T_spin

Select Execute.

Note that the variation of the critical temperature with the DP is almost a perfect straight line. This is consistent with the classical Flory-Huggins theory and with recent experimental and simulation results for isotopic blends (Gehlsen et al. 1992; Deutsch and Binder 1992). If you wish, you may perform a similar calculation using the MSA as the closure to see that the MSA predicts a critical temperature proportional to the square root of the DP. A related consequence is that the values of the critical demixing temperatures predicted by the MSA for UCST systems are far lower than those predicted by the R-MMSA.

As the results of this tutorial demonstrate, the R-MMSA is in principle the preferred closure. However, difficulties in solving the R-MMSA equations may make the MSA the preferred closure in many cases. For melts, the MSA should be adequate. For blends in which a trend is of greater interest than the precise value of Tc, the MSA also may be useful.

Finally, recall that in the R-MMSA, the hard core system is used as a reference system, which is solved with the MSA. Using current methods of solution, it is not possible to use the R-MMSA when the reference system is phase-separated. Thus, computing LCST phase diagrams using the method of Lesson 5 given online through Pilot is not now possible with the R-MMSA.

Use the Session/Quit command to exit Insight II.




Last updated April 20, 1998 at 09:59AM PDT.
Copyright © 1997, Accelrys Inc. All rights reserved.