PRISM Theory Study of Amphiphilic Block Copolymer Solutions with

Sep 7, 2017 - We present a comparison of Polymer Reference Interaction Site Model (PRISM) theory and molecular dynamics (MD) simulations for studying ...
0 downloads 10 Views 4MB Size
Article pubs.acs.org/Macromolecules

PRISM Theory Study of Amphiphilic Block Copolymer Solutions with Varying Copolymer Sequence and Composition Ivan Lyubimov, Daniel J. Beltran-Villegas, and Arthi Jayaraman* Department of Chemical & Biomolecular Engineering, Colburn Laboratory, University of Delaware, 150 Academy Street, Newark, Delaware 19716, United States S Supporting Information *

ABSTRACT: We present a comparison of Polymer Reference Interaction Site Model (PRISM) theory and molecular dynamics (MD) simulations for studying amphiphilic block copolymers in solution. We use a generic coarse-grained model to represent amphiphilic A−B block copolymers in implicit solvent with the solvophobicity of the B segments captured using effective B−B pairwise attraction modeled using the Lennard-Jones potential. We study the assembly of the amphiphilic A−B block copolymer as a function of solvophobicity for varying copolymer sequences (diblock and triblock) and composition (solvophobic-rich or solvophilic-rich). The PRISM theory equation along with the atomic Percus−Yevick closure is solved to obtain the intermolecular pair correlations in real space, g(r), and structure factors in Fourier space, S(k), for block copolymer solutions at increasing values of solvophobicity. The real-space intermolecular pair correlation functions and structure factors from PRISM theory and from MD simulations are compared directly. We find excellent quantitative agreement in g(r) between PRISM predictions and MD simulations at low solvophobicities where the block copolymer solution is in a disordered state. PRISM theory captures the concentration fluctuations at low solvophobicities well but fails to converge to a numerical solution at higher solvophobicities where we see evolution of ordered structures in MD simulations. Despite this drawback, PRISM theory is a valuable tool as the low solvophobicity results from PRISM predict many of the thermodynamic and structural signatures of the block copolymer solutions at higher solvophobicities. For example, we find that the inverse microphase peak 1/S(k*) at low solvophobicities obtained from PRISM theory, which when extrapolated to zero quantifies the spinodal transition solvophobicity value, is in good quantitative agreement with MD simulations. Additionally, in those systems where MD simulations predict aggregation of the micelles/clusters at high solvophobicity, the structure factors from PRISM theory at low solvophobicity also present an increasing value in the zero-wave vector structure, indicating a tendency toward macrophase separation at higher solvophobicity. These results show the capability of PRISM theory to predict assembly over a wide range of design parameters of copolymers and guide the use of computationally intensive molecular simulations.

1. INTRODUCTION Polymer Reference Interaction Site Model (PRISM) theory is a liquid state theory developed by Schweizer and Curro1,2 as an extension of RISM theory by Chandler and Andersen.3 PRISM theory has been applied to a broad class of soft macromolecular materials, including homopolymer melts, 2,4,5 polymer blends,6−10 polymer nanocomposites (see references in review article11), and melts and solutions of block copolymers,12−16 providing a platform to predict structure and thermodynamics and interpret experimental measurements on analogous systems.5,15,17−19 PRISM theory has also been directly compared with Monte Carlo (MC) and molecular dynamics (MD) computer simulations to test the ability of PRISM theory to quantitatively or qualitatively predict structural and thermodynamic results seen in molecular simulations.9,18−25 Such direct comparisons between PRISM and simulations have been made for polymer melts, 18 ,20 ,23 ,2 6,2 7 polymer blends,9,19,22,24,25 polymer nanocomposites,11,28,29 and block copolymers.21,30 The focus of this paper is to evaluate the capabilities and limitations of PRISM theory specifically for studying assembly © XXXX American Chemical Society

of amphiphilic block copolymers in solution conditions as a function of varying solvophobicity of one of the blocks. Amphiphilic block copolymers are macromolecules with two or more blocks with different affinities to the solvent/solvent mixture. Experimentally, solvophobicity of one of the blocks in the copolymer can be varied either by changing the temperature of the solution and altering the solvent−block interactions or by changing the composition of two or more solvents in the solution and, in turn, affecting the net solvency of one of the blocks. Past experiments, theory, and simulations have shown that the amphiphilicity or difference in solvent affinity of the blocks drives block copolymer chains to assemble and form ordered/aggregated structures (e.g., micelles, gels, layered aggregates).31−38 In the simplest case of symmetric diblock copolymers the solvophobic blocks tend to aggregate into a spherical core and are shielded from the solvent by solvophilic block coronas. The vast available design space Received: July 4, 2017 Revised: August 23, 2017

A

DOI: 10.1021/acs.macromol.7b01419 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

there is a good agreement between the structures seen in experiments60 and in simulations.46,61 In the case of triblock sequences with solvophilic center block, both simulations and experiments show that the assembled micelles can further aggregate due to bridges formed by the triblock sequence.41,42,62−64 In summary, there is significant literature suggesting good qualitative, and in many cases quantitative agreement, between experiments and molecular simulations. This gives credibility to the use of simulations in this paper for evaluation of PRISM theory. Even though molecular simulation techniques are widely used for block copolymers, limitations in terms of the computational expense (i.e., hours on many nodes of a high performance supercomputer) required to capture the experimentally relevant length and time scales have led to studies of block copolymers using other theoretical approaches like field theoretic polymer simulations, 65 self-consistent field theory,66−70 classical density functional theory (DFT),71,72 and statistical associating fluid theory (SAFT)-based methods.73−76 Each of these approaches has its advantages and limitations in terms of their applicability for studying block copolymer solutions; we direct the reader to the above references for more details of these methods as a review of all of these methods is outside the scope of this paper. In this paper, our primary goal is to evaluate how well PRISM theory predicts the trends in structure and thermodynamics in dilute solutions of A−B amphiphilic block copolymers. We study linear A−B amphiphilic block copolymers with varying copolymer compositions, sequence, and solvophobicity of the solvophobic B block and conduct a quantitative comparison between PRISM and MD simulations for these systems. We calculate the real and reciprocal space pair correlation functions and structure factors, at system design parameters and conditions where PRISM theory provides a numerical solution. We then compare the pair correlation functions and the value of the critical solvophobicity required for block copolymer assembly computed from PRISM theory with that from MD simulations. We find that the results from PRISM and MD show excellent quantitative agreement in the pair correlation functions at the low solvophobicities where PRISM theory converges to results and the system is mostly disordered. Additionally, the features of the structure factor at finite wave vector k* and k → 0 obtained from PRISM theory at low solvophobicities (disordered states) predict quantitatively the phase behavior seen in MD simulations at higher solvophobicities (ordered states). Lastly, the critical solvophobicity needed for assembly of copolymers, approximated to be the spinodal instability, is in good agreement with the critical solvophobicity values from MD simulations. These qualitative and quantitative agreements between PRISM and MD simulation results show that PRISM theory can be used as an effective tool to guide MD simulations and experiments in the study of block copolymer assemblies for various block copolymer sequences and compositions and potentially nonlinear architectures. Additionally, since PRISM calculations are much faster than full MD simulations, it allows the exploration of a much larger parameter space, namely different block copolymer architectures, sequences, compositions, and solvophobicities, with and without nanoscale additives to identify systems and conditions of interest for further study. The rest of the paper is organized as follows. In the Method section, we describe the model used to represent the linear amphiphilic block copolymers in PRISM theory and MD

comprised of not only the chemistry of monomers but also sequence, composition, and architecture leads to a broad range of structures formed by amphiphilic block copolymers in solution, namely spherical39,40 and worm-like micelles39,40 as well as networks/gels.41,42 These assembled structures have found use in many different applications, such as biological and nonbiological cargo delivery,35,43−45 templating inorganic nanoparticle synthesis,46,47 stimulus-responsive gels,42 and emulsifiers/stabilizers.40 In the past, PRISM theory has been applied to study block copolymers of arbitrary composition, sequence, and architecture in solutions as well as in melts.12−16,21,48,49 Early studies by David and Schweizer focused on symmetric and asymmetric block copolymers12−14 in melt-like concentrations using the Gaussian thread model (also called thread PRISM).50 They quantified how chain length, density, temperature, and composition impacted the effective χ parameter and peak scattering intensity. They also compared results from PRISM theory with reference molecular mean spherical approximation (RMMSA) closure to expected behavior based on Leibler’s mean-field analysis51 as well as results from field theoretic treatment of mean-field theory with fluctuation corrections.52 Going beyond the melt-like conditions, Guenza and Schweizer focused on block copolymers in dilute to concentrated solutions, using PRISM theory with Gaussian thread model with molecular closures.15,16 In the case of diblock copolymer solutions under neutral solvent conditions, they predicted local and microdomain scale concentration fluctuations and their impact on the link between degree of polymerization, temperature, and polymer concentration at the order−disorder transition (ODT). They showed that in the semidilute regime PRISM theory predictions agree with blob scaling and fluctuation-corrected field theoretic analyses, while strong disagreements occurred in the concentrated solution and melt regime.16 Nevertheless, at concentrated conditions, they provided predictions for apparent scaling laws with effective exponents being dependent on solvent quality. Going beyond the thread model, PRISM theory has also been used with freely joined chain (FJC) model which is a more realistic model than the Gaussian thread model but needs to be solved numerically. For example, PRISM theory with FJC model has also been used to study block copolymers and compared against Monte Carlo simulations.21 The qualitative dependence of the inverse peak structure factor on temperature and peak wave vector dependence on temperature have been compared. Additionally, PRISM theory studies using molecular closures that were previously tested on polymer blends53 have been used for block copolymer systems to predict polymer length dependence of inverse peak structure factor.52 While the past PRISM theory-based studies of diblock copolymer systems are valuable, to the best of our knowledge, PRISM theory has not been applied to predict structure and thermodynamics during micellization/aggregation in amphiphilic block copolymers for varying copolymer sequence and composition and solvophobicity, which is the focus of this paper. We note that there have been many simulations studies that have captured the solution behavior of block copolymers for varying sequences, compositions, flexibilities, etc. For example, Monte Carlo (MC),54 molecular dynamics (MD),55 and dissipative particle dynamics (DPD) simulations56,57 reproduce experimentally39,58,59 observed micellar morphologies obtained for varying diblock copolymer composition. In the case of triblock sequences with solvophobic center block, B

DOI: 10.1021/acs.macromol.7b01419 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules simulations, the various design parameters, and conditions explored, namely block copolymer compositions and sequences, solvophobicities, etc., and a brief overview of PRISM theory and MD simulations along with the data analyses. In the Results and Discussion section, we compare intermolecular correlations between PRISM and MD simulations in real and reciprocal space and discuss thermodynamics and PRISM predictions for phase transition conditions as compared with data from MD simulations. Finally, in the Conclusion section we provide a summary of our results, implications of our findings, and our vision for future work.

Table 1. Range of A−B Copolymer Systems Studied in This Paper composition

description

notation

sequence

symmetric

diblock triblock (B-center) triblock (A-center)

case 1 case 2 case 2i

A12-b-B12 A6-b-B12-b-A6 B6-b-A12-b-B6

diblock diblock (inverse) triblock triblock (inverse)

case case case case

A6-b-B18 B6-b-A18 A3-b-B18-b-A3 B3-b-A18-b-B3

asymmetric

2. METHOD 2.1. Model. We model the amphiphilic A−B copolymers of interest using a coarse-grained (CG) model where each CG bead of diameter d represents a Kuhn segment of A and B chemistry. While PRISM theory (described in section 2.3) uses a CG model with a fixed bond length between bonded neighboring beads, the molecular dynamics simulation method (described in section 2.4) uses a bead−spring model77,78 with an equilibrium distance of 1d and a spring constant of 50ε/d2, where ε is the energy unit in the simulation. The effect of the solvent is modeled implicitly through the interaction potential between nonbonded pairs of A and B CG beads. The A−A and A−B nonbonded pairwise interactions are modeled using the purely repulsive Weeks−Chandler−Andersen (WCA) potential79 with a cutoff distance of 21/6d, and B−B nonbonded pairwise interactions are modeled using Lennard-Jones (LJ) potential. While the B−B LJ potential cutoff distance is maintained constant at 2.5d, the attraction strength of the B−B LJ potential, εBB, is varied to capture the increasing solvophobicity of B segments or monomers in the implicit solvent. The implicit solvent model where the solvophobicity is captured through effective attraction between polymer segments is equivalent, in a thermodynamic sense, to explicit model descriptions of polymer solutions where the solvent quality is captured by pairwise interactions between polymer segments and solvent molecules. Molecular simulation studies that compare implicit and explicit solvent descriptions of solutions of homopolymers80 and copolymers81 have shown that implicit solvent models can be mapped to mimic the polymer conformations and structure seen in the explicit solvent models due to the solvent quality. The comparison of relaxation dynamics has shown that implicit solvent models are more prone to be stuck in metastable states than explicit solvent models. In our work we seek to compare equilibrium (or close to equilibrium) correlations between PRISM theory and MD simulations of block copolymer micelles rather than assess the dynamic aspects of the micellization process. Through the protocol we use to slowly increase solvophobicity, we ensure the system is as close to equilibrium as possible (as described in section 2.4), so the differences in relaxation dynamics between implicit and explicit solvent models are not a concern. 2.2. Systems Studied. In this paper, we focus on linear copolymer chains in an implicit solvent at a total volume fraction of η = 0.1, representing a semidilute solution. The A−B block copolymers are at constant chain length equal to 24 CG beads, with varying sequences (diblock or triblock) and compositions (symmetric or asymmetric, either A-rich or Brich). In Table 1, we list all the systems we study in this paper. Within the symmetric copolymers, with equal number of type A

3 3i 4 4i

and type B beads, we study two sequences: diblock (case 1) and triblock copolymers with either solvophobic (B-type) or solvophilic (A-type) center block, labeled as case 2 and case 2i (or short for case 2 inverse), respectively. Within the asymmetric copolymers, we limit ourselves to the four cases summarized in Table 1. 2.3. PRISM Theory. PRISM theory2,12,48,82 is an integral equation theory described as H(r ) =

∫ dr′ ∫ dr″ Ω(|r − r′|)C(|r′ − r″|)[Ω(r″) + H(r″)] (1)

or in Fourier space expressed as Ĥ (k) = Ω̂(k)Ĉ(k)[Ω̂(k) + Ĥ (k)]

(2)

where Ĥ (k) is a N × N matrix representing the intermolecular pair correlations, Ω̂(k) is an N × N matrix representing the intramolecular pair correlations, and Ĉ (k) is an N × N matrix representing the direct pair correlations for a system with N types of sites. In our study of A−B block copolymers we have two types of sites, A and B, arranged in varying sequences (diblock versus triblock) and compositions (majority A to minority A) leading to 2 × 2 matrices in eq 2. The components of Ω̂(k), ωαγ(k), capture the shape of the polymer chain and is a known input to the PRISM equation. Early implementations of PRISM theory used analytic expressions for ωαγ(k) (e.g., ideal thread chain12 or freely jointed chain models8,83). Alternately, ωαγ(k) can be obtained by sampling configurations of the molecule in the conditions of interest in molecular simulations.29,48 In this paper, we use trajectories from two cases of molecular dynamics (MD) simulations to sample A−A, A−B, and B−B intramolecular pair correlations, and calculate ωαγ(k) as ωαγ (k) =

1 Nα + Nγ





∑∑ i=1 j=1

sin(krij) krij

(3)

where Nα and Nγ are the number of sites per chain of types α and γ, rij is the distance between sites i and j, and angle brackets represent statistical averaging. In one case, we use the MD trajectory of a single A−B copolymer chain in implicit solvent (denoted as single-chain MD, scMD) as the source for ωαγ(k). This scMD ωαγ(k) input neglects the effects from other chains at the desired concentration on the intramolecular correlation. In a second case, we calculate ωαγ(k) from the trajectories obtained from computationally intensive MD simulations (denoted as full-fledged MD, ffMD, purely to contrast against scMD), with the exact number of polymer chains present in the simulation box needed to achieve solution volume fraction of η = 0.1. We compare the results from PRISM using scMD and C

DOI: 10.1021/acs.macromol.7b01419 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules ffMD ωαγ(k) as inputs to see if one is better than the other in reproducing analogous results from a MD simulation. We also note that there have been studies where the ωαγ(k) is obtained in a self-consistent framework where the solvation potential calculated from PRISM theory is used in the molecular simulation of a single molecule to iterate from an initial guess of the ωαγ(k) to the ωαγ(k) of the molecule in the conditions of interest.18,22,28,84−86 We choose not to adopt the self-consistent approach in this paper as we see excellent agreement in the disordered state structure between PRISM theory results and simulations using the scMD and ffMD ωαγ(k). While this is the case for the systems and conditions we present in this work, it may not necessarily be so in other systems and conditions where the influence of other chains/molecules in the system on the ωαγ(k) of the molecule warrants a self-consistent approach. In such cases, either an analytical expression of ωαγ(k) (e.g., ideal conformation) or ωαγ(k) obtained using scMD will need to be iterated through the self-consistent PRISM-MD approach to obtain correct results. Since the desired outputs, Ĥ (k) and Ĉ (k), are two unknown matrices, to solve the PRISM equation, we need a closure relationship. The closure relates the direct pair correlation functions Cαγ(r), intermolecular pair correlation hαγ(r), defined as hαγ(r) = gαγ(r) − 1, and the pairwise interaction potential, Uαγ(r). Here we choose to use the atomic Percus−Yevick (PY) closure87 for A−A, A−B, and B−B pairs, which in real space can be written as

pairwise structure factor, S(k), is calculated using the intramolecular and intermolecular pair correlation matrices as follows: S(k) = Ω̂(k) + Ĥ (k)

The diagonal components of the static structure factor matrix can also be calculated as Sαα(k) = 1 + 4πρα

r≤d

∫0



tot r dr(gαα (r ) − 1)

sin(kr ) k

(6)

gtot αα(r)

where is the total radial distribution function (having both inter- and intramolecular pair correlations). 2.4. Molecular Dynamics (MD) Simulations. We conduct NVT molecular dynamics simulations using the LAMMPS package.96 We create an initial configuration of the system of copolymer chains with each chain in a rod-like conformation, on a cubic lattice at a low η, high temperature, and B−B nonbonded LJ interaction strength, εBB = 0.1. The total number of copolymer CG beads in each simulation is 14 400, corresponding to 600 copolymer chains of length 24 beads each. Periodic boundary conditions are imposed on all directions. We compress this initial simulation cubic box size isotropically to a size that results in the desired volume fraction η of 0.1; we also relax the chains from their unphysical rod-like initial configuration. We then lower the T* over 600 000 time steps, where each time step corresponds to Δt = 0.005, from an initial value of T* = 4 to T* = 1. Since the system is fluid-like at this low εBB = 0.1, we expect that the annealing schedule will not affect the eventual results. After an additional 600 000 time steps at T* = 1 we begin annealing the system through a carefully selected stagewise increase in εBB. In each stage, εBB is kept constant, and the simulation is run for 106 time steps, and then to go to the next stage the εBB is increased by ΔεBB = 0.009. This staged increase of εBB mimics the exploration of copolymer chain equilibrium assembly in solution as a function of solvophobicity captured by εBB. In each stage, the simulation box configurations are saved every 100 000 time steps. Our choice of annealing schedule quantified by the stage length, i.e., the number of time steps at each stage before increasing εBB and the stage height, i.e., ΔεBB, is chosen to ensure we have sufficient time for equilibration in each stage. Our choice of sampling frequency ensures statistical independence between sampled configurations. To compare to the outputs from PRISM theory, in MD simulations we calculate the real-space pair correlations gαγ(r) using standard techniques97 as well as Fourier space structure factors.29 We also identify micellar clusters by grouping copolymer chains that are part of the same cluster. We define any two chains to be part of a cluster if they satisfy the “connectivity criterion”. The “connectivity criterion” is defined as follows depending on the copolymer sequence and composition: Connectivity criterion for cases 1, 2, 3, 3i, and 4: We identify two copolymer chains as neighbors if the distance between their respective solvophobic B block centers of mass is below a coordination threshold. The coordination threshold is the first valley of the solvophobic B block center of mass radial distribution, gBB(r), at high εBB = 1.0; this is found to be equal to 5d. Two chains are part of the same cluster and satisfy the connectivity criterion, if these two chains are neighbors and if each of these chains has at least a minimum number of neighbors. The minimum number of neighbors is identified by

Cαγ(r ) = (1 − e βUαγ(r))gαγ (r ), r > d gαγ (r ) = 0,

(5)

(4)

Past work has looked into the validity of atomic closures and molecular closures in systems with two or more distinct sites.27,53,88,89 In general, in macromolecular systems with strong attractive interactions (e.g., strongly segregating block copolymers or blends), molecular closures are expected to perform better than atomic closures. For example, in the case of phase separating polymer blends, it was found that molecular closures predict the chain length dependence of the critical temperature53,88 better than atomic closures do. In general, the conclusion from these studies was that atomic closures are better suited for small molecules than macromolecules, and in the case of macromolecules, the atomic closures are better suited for systems with athermal interactions than strongly attractive interactions (e.g., strongly segregating block copolymers). Recent work using PRISM theory with atomic closures to study polymer nanocomposites with athermal/weakly attractive interactions has shown good agreement with results from MD simulations and experiments.90−92 Given the recent success with atomic closures studying systems with weak attractions and the reduced complexity in the numerical convergence of atomic closures compared to molecular closures, we choose to use atomic closures in this work to test the ability of PRISM theory to predict structure and thermodynamics of micellization in copolymer solutions. We numerically solve the coupled nonlinear integral equations (2) and (4) using the KINSOL algorithm.93 The KINSOL algorithm produces numerical solutions faster and for a larger design space than the Picard technique used in many previous PRISM papers.84,90,94,95 Upon solving, the two outputs we obtain are intermolecular pair correlation functions, gαγ(r), and structure factors, Sαγ(k). The 2 × 2 matrix for the D

DOI: 10.1021/acs.macromol.7b01419 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

Figure 1. Intermolecular correlation function, gαγ(r), from PRISM theory (lines) using the two different intramolecular structure factors (scMD and ffMD) and from MD simulations (symbols) for symmetric blends (green circles), case 1 (blue squares), case 2 (red diamonds), and case 2i (orange triangles). In each plot, to clearly show the g(r) curves one above the other we shift each subsequent curve up by 0.6. Left, center, and right panels correspond to AA, AB, and BB components of gαγ(r), respectively. The lines on the top row (scMD) correspond to PRISM results where single chain MD is used to obtain intramolecular distribution function, ωαγ(k). The lines on the bottom row (ffMD) correspond to PRISM results where ωαγ(k) is from MD simulations with systems at packing fraction of η = 0.1. The results in Figure 1 are at a constant solvophobicity parameter εBB = 0.1. n

neighbor distribution functions within the coordination threshold and is used as an additional criterion to distinguish between copolymer chains in fluid-like environments as opposed to copolymer chains in clustered state. Connectivity criterion for cases 2i and 4i described in Table 1: The topological separation of the solvophobic B blocks in these two cases causes the beads of these two solvophobic end blocks to either belong to the same cluster or micelle64,98 or belong to two different micelles with the solvophilic A block serving as the bridge between the micelles.64,98 We identify micellar clusters formed by cases 2i and 4i chains by splitting each chain in two equal parts and proceed to use similar connectivity criteria as for case 1 using the centers of mass of each B block separately. We consider two B block ends to be neighbors if the distance between their centers of mass is below a coordination threshold of 5d; furthermore, the B block ends are “connected” if each solvophobic end has a minimum of 13 end neighbors. We follow the same rationale to determine connectivity criteria parameters for cases 2i and 4i as done for cases 1, 2, 3, 3i, and 4. Additional details regarding the determination of connectivity criteria can be found in the Supporting Information section S.I and Figures S1 and S2. After pairs of chains (or pairs of B block ends) are identified as part of the same micellar cluster, we use a “friends-of-friends” algorithm99 to identify all clusters. Once all the clusters are identified, the number of clusters (of sizes greater than 1) denoted as n and the size of each cluster denoted as Nclus,i for ith cluster is recorded. Using the n number of clusters and simulation box volume Vtot, we calculate the cluster density/ concentration as n ρclust = Vtot (7)

⟨Nclus⟩ =

∑i = 1 Nclus, i n

(8)

3. RESULTS AND DISCUSSION 3.1. Pair Correlations in Real Space. First, we compare the gαγ(r) for the symmetric cases in Table 1 from PRISM theory using either scMD or ffMD intramolecular structure factors, as described in the Method section. For all cases, we also present the gαγ(r) obtained directly from MD simulations. In Figure 1, for all symmetric copolymer cases and all pairs of gαγ(r), we see excellent quantitative agreement between the PRISM theory results using scMD and ffMD intramolecular structure factors. It is not surprising that the intermolecular pair correlation functions obtained using PRISM theory with scMD and ffMD intramolecular structure factors are in agreement as the corresponding input scMD and ffMD ωαγ(k) are also in agreement (see Supporting Information Figures S3 and S4). The scMD and ffMD ωαγ(k) for the systems presented in Figure 1 at the solution concentration η = 0.1 are almost indistinguishable at small wave vectors and have negligible differences at large wave vectors. The results in Figure 1 are at a constant solvophobicity parameter εBB = 0.1, which is the highest value of solvophobicity parameter εBB for which PRISM theory converges to numerical solution for the blends case. Even though for the other copolymer cases in Table 1, PRISM converges to a numerical solution at higher εBB values to enable comparison between all cases we choose to present these four systems at the solvophobicity parameter εBB = 0.1. A comparison of gαγ(r) at higher εBB for these symmetric cases is presented in Figure S5. We also see good agreement between PRISM results from scMD and ffMD for all symmetric copolymers at the higher εBB where the system still maintains a disordered state.

Using the individual cluster sizes and number of clusters, we also calculate the average cluster size as E

DOI: 10.1021/acs.macromol.7b01419 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

Figure 2. Intermolecular correlation function, gαγ(r), from PRISM theory (lines) using the two different intramolecular structure factors (scMD and ffMD) and from MD simulations (symbols) for case 3 (red circles), case 3i (orange squares), case 4 (blue diamonds), and case 4i (cyan triangles). In each plot, to clearly show the g(r) curves one above the other we shift each subsequent curve up by 0.6. Left, center, and right panels correspond to AA, AB, and BB components of gαγ(r), respectively. The lines on the top row (scMD) correspond to PRISM results where single chain MD is used to obtain intramolecular distribution function, ωαγ(k). The lines on the bottom row (ffMD) correspond to PRISM results where ωαγ(k) is from MD simulations with systems at packing fraction of η = 0.1. The results in Figure 2 are at a constant solvophobicity parameter εBB = 0.1.

is an expected result as the A-block is modeled as athermal, and A−A and A−B pair correlations are not directly affected by change of εBB. There is a noticeably stronger B−B correlation in all cases for higher εBB with the appearance of the first nearestneighbor shell in gBB(r). As seen for εBB = 0.1, the B−B correlations at εBB = 0.27 in case 2 are smaller than for case 2i due to the sequence-originated “screening” effect as discussed in the previous paragraph. In Figure 2, the gαγ(r) for asymmetric cases is presented. As in the symmetric cases, for the asymmetric cases (Figure 2) we also see a quantitative agreement between the results from PRISM using either inputs (scMD or ffMD) and MD simulations at εBB = 0.1. Comparing the B−B pair correlations, the inverse cases (cases 3i and 4i) have smaller correlation holes than cases 3 and 4. This is because cases 3 and 4 have larger solvophobic B block (majority composition is that of B monomers) than their inverse counterparts, leading to larger intramolecular B−B contacts and, in turn, a larger correlation hole in the intermolecular correlation. Additionally, in case 4 the intermolecular B−B interactions are “screened” by the presence of A-blocks at the ends of the chain. It can be also seen from the gAA(r) that for case 4 the range of correlation hole is shorter than for other three cases present in the figure. This is because unlike the other three cases, the case 4 sequence naturally favors intermolecular A−A interaction the most by having A-blocks at both ends of the triblock copolymer chains. So far, in Figures 1 and 2 we see no differences between the results from the two different ωαγ(k) (from scMD and from ffMD). Since there is no noticeable effect of including many chain interactions into intramolecular distribution for the systems and conditions (disordered block copolymer solutions) studied here, in the remainder of this paper, we use scMD as the way to obtain input ωαγ(k) for PRISM results. This is important as it allows us to capitalize on the fact that the scMDbased PRISM calculation takes only a small fraction of the computational time required for ffMD-based PRISM. More

In Figure 1, comparison of the results from PRISM (scMD or ffMD) with that from MD simulations shows an overall good quantitative agreement at solvophobicity parameter εBB = 0.1. The minor differences in gαγ(r) between PRISM and MD at low r is not surprising as they have also been seen in previous studies on polymer systems aimed particularly at comparing PRISM theory and simulations for homopolymer melts and polymer blends.5,9 The sharp edges/cusps in the PRISM gαγ(r) curves at distances around r = 1 are due to discontinuous closure condition, eq 3, and in principle can be smoothed upon usage of other closure relations. If we compare the gαγ(r) (intermolecular pair correlation) for all the cases at this solvophobicity of εBB = 0.1, we see a correlation hole at short distances commensurate with the size of the blocks in the same chain (in the case of copolymers) or size of chains (in the case of blends). There are some differences between the cases, such as a stronger B−B pair correlation in the case of symmetric blends than symmetric copolymers. This is an expected result as the covalently linked A-block in the case of copolymers constrains B−B intermolecular contacts; this constraint is absent in the blends. Another difference is the smaller correlation hole in A−A pair correlations for case 2 than case 2i and the opposite trend for B−B correlations. In case 2, the middle solvophobic B-block is surrounded by two solvophilic A-blocks, which leads to “screening” of some B−B intermolecular interactions. In case 2i, the solvophobic B-blocks located at the chain ends favors unrestrained B−B intermolecular contacts. PRISM theory clearly captures these molecular level packing trends quantitatively correctly for these symmetric block copolymer solutions. The effect of increased solvophobicity on intermolecular pair correlations in symmetric block copolymers is shown in Figure S5. For all symmetric block copolymer cases, the A−A and A− B pair correlations are unaffected by the increase of solvophobicity parameter from εBB = 0.1 to εBB = 0.27. This F

DOI: 10.1021/acs.macromol.7b01419 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

macrophase separation. This emergence of a microphase peak even at low solvophobicities, where the solution is effectively disordered, is analogous to the “fluctuating microstructures” mentioned in previous PRISM studies by Guenza and Schweizer” and was also seen for compositionally asymmetric diblock copolymers far from the order−disorder transition,16 in agreement with experimental observations by Lodge and coworkers.100 Similar observations have been discussed in other past work on symmetric and asymmetric diblock copolymer systems.39,40,101,102 For case 2, A6-b-B12-b-A6, in Figures 3b and 3e we see a more subdued k* peak growth in SBB(k) as compared to case 1 for the same solvophobicity values. This is because the A-blocks at either end of the solvophobic B block sterically hinder the B−B intermolecular contacts and reduce concentration fluctuations. For case 2i, B6-b-A12-b-B6, in Figures 3c and 3f we see behavior not seen in case 1 or case 2. The SAA(k → 0) shows an upturn and a tendency to diverge as εBB increases, and SBB(k) shows increasing microphase peak as well as increasing SBB(k → 0). This suggests that as solvophobicity is increased, the disordered solutions have concentration fluctuations at multiple length scales, such that at higher solvophobicity (beyond the values shown in Figure 3) we can expect the system with case 2i polymers to show microphase order and expect those microphase-separated structures (e.g., micelles) to aggregate and macrophase separate. Past work focused on such sequences have shown that the two solvophobic end blocks at high solvophobicity could bridge the microphase-separated domains and cause gelation/macrophase separation.41,63 The results in Figure 3 are for those εBB values at which PRISM theory converged to a numerical solution; at higher εBB PRISM theory did not converge to a solution. In contrast to PRISM theory, MD simulations have no such convergence issues and are used to simulate these systems at larger εBB values. We present the simulation snapshots from the MD simulations in Figure 4 that confirm what we predict at high εBB based on the low εBB PRISM results. Figure 4 presents the simulation snapshots at low εBB where PRISM theory converged (first column), an intermediate εBB (middle column), and a high εBB (right column); PRISM theory did not converge to a numerical solution at both the intermediate and high εBB. The specific values of intermediate and high εBB are chosen based on some key structural and thermodynamic transition points we see in MD simulations as described in detail in section 3.3. For each copolymer case, the intermediate εBB represents the highest value of εBB where we observe an isotropic structure in that copolymer system. Beyond that the intermediate εBB we see ordering in the solution and in some cases anisotropic structures forming. For cases 1 and 2, as predicted by the PRISM structure factors in Figure 3, we see in the snapshots microphase-separated micellar structure and no macrophase separation as εBB increases. For case 2i, we see concentration fluctuations and microphase-separated structures which then aggregate at higher εBB as predicted by PRISM theory structure factors at lower εBB. In addition to qualitatively predicted tendency of these copolymer systems to microphase and/or macrophase separate at higher εBB, PRISM theory can also provide some quantitative predictions at higher εBB confirmed with MD simulations snapshots and analysis. First, the k* denotes the length scale of the emerging concentration fluctuations and in real space corresponds to r = 2π/k*. Based on the k* being the smallest for case 1, we expect the aggregate/micelle structures for case 1

importantly, the computational time needed to complete an scMD-based PRISM calculation is significantly smaller than simulating an analogous system using our MD simulation protocol. For example, for case 1, the symmetric diblock copolymer solution, the wall time (equilibration and production) to complete an MD simulation using the protocol described in the Method section over the range of solvophobicities εBB = 0.1−1.0 (going from disordered fluid to ordered micelles) is 34 h on 4 CPU cores in the University of Delaware HPC Farber cluster composed of Intel(R) Xeon(R) CPU E5-2670 v2 @ 2.50 GHz processors. For the PRISM calculations over the range of solvophobicities εBB = 0.1−0.4 (the highest value where numerical convergence is possible), the time needed for the scMD simulation and calculation of ω(k) and numerical solution of PRISM equations together is ∼4.5 min on all 4 cores of the Intel Xeon E3-1225 v3 Quad Core, 3.20 GHz Turbo processor. Comparing the time for MD simulation and PRISM calculation for the same range of solvophobicities εBB = 0.1−0.4, we estimate that the PRISM calculation is ∼78 (= 5.88 h/4.5 min) times faster than the corresponding MD simulation. 3.2. Structure Factors. In section 3.1 we discuss real-space intermolecular pair correlation for the various sequences and compositions in Table 1. While these real-space intermolecular pair correlation functions are useful in understanding localized arrangement of monomers, to better understand how concentration fluctuations emerge and lead the system toward assembly and phase separation, in this section we show the A− A and B−B structure factors for all the cases. In Figure 3, for case 1, A12-b-B12, as εBB increases the SAA(k) shows little variation, while SBB(k) exhibits increasing value at the small-angle wave vector k* (marked with the arrow in Figure 3d), indicating concentration fluctuations leading to an emergence of a microphase peak and the tendency to microphase separate. We do not see any increase in the structure factor as k → 0, confirming lack of tendency for

Figure 3. Static (a−c) A−A and (d−f) B−B structure factors from PRISM theory with ωαγ(k) from scMD are shown for symmetric cases. The results are for case 1 (A12-b-B12), case 2 (A6-b-B12-b-A6), and case 2i (B6-b-A12-b-B6). We present only the structure over a small range of wave vectors in this figure, even though (as seen in the inset in the top right panel) the structure factor output from PRISM covers a broader range of wave vectors. The values of solvophobicity, εBB, are shown in the legend. G

DOI: 10.1021/acs.macromol.7b01419 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

case 2, and case 2iat one specific solvophobicity εBB = 0.25 to demonstrate the quantitative agreement in k* values between PRISM theory and MD simulations. The theoretical lines are almost indistinguishable from simulation data for k > 1.2, whereas deviations occur at smaller wave vectors. These deviations are more clearly seen in Figure S6 and originate in the inability of MD simulation with a finite simulation box to capture accurately large scale (small k) structure. This also shows the value of PRISM in evaluating the accuracy of simulations at large scales and signaling when the box sizes are likely to be too small or incommensurate with the ordering periodicity. We note that besides some differences between absolute values of SAA(k) and SBB(k) at small k from PRISM and MD, we see a good agreement in the location of k* values for the peaks. In Figure 6, we present A−A and B−B structure factors with increasing solvophobicity for cases 3, 3i, 4, and 4i. In all cases,

Figure 4. Snapshots from MD simulations for symmetric cases, where A-type beads are shown in red and B-type beads are green. The snapshots are for (a) case 1 (A12-b-B12), (b) case 2 (A6-b-B12-b-A6), and (c) case 2i (B6-b-A12-b-B6) at low εBB, corresponding to a liquidlike structure (left), at intermediate εBB (center), and at high εBB, corresponding to microphase-separated and/or aggregated structures (right).

to be larger than those compared to cases 2 and 2i. Figure 4 simulation snapshots confirm this behavior, especially at the highest solvophobicity. Figure 5 and zoomed-in images in Figure S6 present the direct comparison of S(k) from PRISM theory and from MD simulation for all three casescase 1,

Figure 6. Static (a−d) A−A and (e−h) B−B structure factors from scMD-based PRISM theory for asymmetric cases, namely case 3 (A6-bB18), case 3i (B6-b-A18), case 4 (A3-b-B18-b-A3), and case 4i (B3-b-A18b-B3). The values of solvophobicity, εBB, are shown in the legend.

we see the tendency for micelle formation with increasing solvophobicity, evident from the growing microphase peak. Among all cases in Figure 6, the lowest k* value in SBB(k*) is for case 3 (Figure 6e), which corresponds to the largest length scale of concentration fluctuations and prediction of largest micelles at high solvophobicity compared to other cases. Additionally, based on the SBB(k*) values (Figure 6h), case 4i shows the slowest peak growth for the range of εBB considered. This is expected since there are only three B-type beads on each chain end separated by long A-block, and the concentration fluctuations are decreased. Cases 3 and 4 also exhibit growing zero wave vector structure in SBB(k), i.e., tendency for macrophase separation, which is absent in cases 3i and 4i. Figures 7 and 8 show the simulation snapshots for three values of εBB as well as the structure factors from MD simulations and PRISM theory at one particular εBB. The simulation snapshots in Figure 7 confirm the above PRISM-based predictions: (a) the largest micelle structure is observed for case 3, (b) case 4i

Figure 5. Static A−A (a) and B−B (b) structure factors from PRISM theory with ωαγ(k) from scMD (lines) are compared with results from MD simulations (symbols) for symmetric cases: case 1 (A12-b-B12) in blue, case 2 (A6-b-B12-b-A6) in red, and case 2i (B6-b-A12-b-B6) in orange, at εBB = 0.25. The inset in part b shows the structure factor SBB(k) for a larger range of k values than the main plot. H

DOI: 10.1021/acs.macromol.7b01419 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

Figure 8. Static A−A (left) and B−B (right) structure factors from PRISM theory with ωαγ(k) from scMD (lines) and from MD simulations (symbols) for asymmetric cases. The results are for case 3 (A6-b-B18) in red, case 3i (B6-b-A18) in orange, case 4 (A3-b-B18-b-A3) in blue, and case 4i (B3-b-A18-b-B3) in cyan all corresponding to solvophobicity εBB = 0.25.

worm-like micelles that aggregate, while in case 3 spherical micelles form a mixture of spherical and some worm-like micelles. Isotropic structure factors obtained from PRISM theory (and MD simulations shown in Figure 8) do not distinguish the two different macrophase-separated structures. 3.3. Micellization Thermodynamics. Next, we explore the thermodynamic information we can gather from PRISM theory. One can use the microphase order parameter 1/SBB(k*) vs εBB data and extrapolate to zero value to find the εBB where SBB(k*) would diverge.16,51 This has been a way that the microphase spinodal transition has been identified in past work on copolymers and single polymer-tethered nanoparticles16,51,103 and could serve as the transition solvophobicity value or εtrBB in our systems in this work. We do this calculation using 1/SBB(k*) vs εBB data from both PRISM and MD simulations. In Figure 9a, the 1/SBB(k*) calculated from PRISM (filled symbols) is compared with 1/SBB(k*) from MD simulations (open symbols) for the symmetric cases. The range of εBB at which PRISM 1/SBB(k*) is presented is limited to low solvophobicity values because of the lack of numerical PRISM solutions beyond a certain solvophobicity. We can cover a larger range of εBB with MD, and thus, for MD results the 1/SBB(k*) is presented over a larger range of εBB. The dotted lines in Figure 9a show extrapolation of PRISM 1/ SBB(k*) to obtain εtrBB. The values of εtrBB calculated via this extrapolation are tabulated in the first column of Table 2. We also obtain the corresponding εtrBB from MD simulation results in Figure 9a and present that in the second column of Table 2. While the values of the microphase εtrBB from PRISM and MD simulations do not show exact quantitative agreement, the values from the two approaches are close and the qualitative trends in how the spinodal εtrBB changes with the various cases is captured correctly. This is in line with previous results for other materials systems, where analogous analysis was performed to compare PRISM with simulations of single polymer tethered

Figure 7. Snapshots from MD simulations for asymmetric cases, where A-type beads are shown in red and B-type beads are green. The snapshots are for namely case 3 (A6-b-B18), case 3i (B6-b-A18), case 4 (A3-b-B18-b-A3), and case 4i (B3-b-A18-b-B3). For each case, we show three snapshots: at low εBB, corresponding to a liquid-like structure (left), at intermediate εBB (center), and at high εBB, corresponding to microphase-separated or aggregated structures (right).

shows the least clear interfaces between the A and B monomerrich regions in the weakly microphase separated state, and (c) macrophase-separated state for cases 3 and 4 at highest solvophobicity, but not so in cases 3i and 4i. We note that the solvophobic blocks in case 4i are small (only three beads) and topologically separated by a long solvophilic A-block (unlike case 3i which is a diblock copolymer); this makes for a weaker driving force for microphase separation as well as macrophase separation. However, at high solvophobicity, MD simulation of case 4i solution shows bridged micelles; the B monomers form cores that are bridged by A blocks. The increasing SAA(k→0), absence of growing SAA(k*), and increasing SBB(k*) with an absence of growing SBB(k→0) with increasing solvophobicity in the PRISM results (Figures 6d and 6h) seem to predict the simulation observation of microphase-separated B cores with A blocks bridging the cores at higher solvophobicity. So far, we have made the case of PRISM theory’s capability to predict qualitatively and quantitatively the structural arrangements at solvophobicities higher than the solvophobicity that converged to a solution. It is also important to note the structural information that is not obtained from PRISM theory which also justifies the use of MD simulations in conjunction with PRISM theory. For example, in the simulation snapshots in Figure 7, at the highest solvophobicity shown, case 4 forms I

DOI: 10.1021/acs.macromol.7b01419 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

In Figure 10a, the 1/SBB(k*) calculated from PRISM (filled symbols) is compared with 1/SBB(k*) from MD simulations

Figure 9. (a) Microphase order parameter, 1/SBB(k*), obtained from PRISM theory (filled symbols) and from MD simulations (open symbols) and (b) cluster number density in reduced units, ρclust, calculated from MD simulations as a function of increasing εBB for case 1 (A12-b-B12) in blue circles, case 2 (A6-b-B12-b-A6) in red squares, and case 2i (B6-b-A12-b-B6) in orange diamonds. Vertical dashed line in (b) plots correspond to cusp positions of ρclust.

Figure 10. (a) Microphase order parameter, 1/SBB(k*), obtained from PRISM theory (filled symbols) is compared with the 1/SBB(k*) from MD simulations (open symbols) for asymmetric cases. (b) Cluster number density, ρclust, calculated from MD simulations is shown as a function of increasing solvophobicity parameter εBB. The results are for case 3 (A6-b-B18) in red circles, case 3i (B6-b-A18) in orange squares, case 4 (A3-b-B18-b-A3) in blue diamonds, and case 4i (B3-b-A18-b-B3) in cyan triangles. Vertical dashed lines in (b) panels correspond to cusp positions of ρclust.

Table 2. Comparison of the Microphase Transition Solvophobicity, εtrBB, Obtained from Extrapolation of PRISM 1/SBB(k*), Extrapolation of MD 1/SBB(k*), and Cusp of ρclust from Cluster Analysis of MD Trajectorya system

εtrBB (PRISM)

case case case case case case case

0.75 0.88 0.92 0.55 0.79 0.64 1.02

1 2 2i 3 3i 4 4i

± ± ± ± ± ± ±

0.04 0.03 0.10 0.02 0.04 0.02 0.06

εtrBB (MD)

εtrBB (MD, ρclust cusp)

± ± ± ± ± ± ±

0.703 0.946 0.922 0.616 0.886 0.706 1.418

0.65 0.83 0.87 0.58 0.89 0.69 1.28

0.01 0.09 0.07 0.04 0.08 0.08 0.05

(open symbols) for the asymmetric cases. The dotted lines in Figure 10a show extrapolation of PRISM 1/SBB(k*) to obtain εtrBB. The values of εtrBB calculated via this extrapolation is in Table 2, first column. We also obtain the corresponding εtrBB from the MD simulation data in Figure 10a (Table 2, second column) and from the cusp of ρclust in Figure 10b (Table 2, third column). The values from the three columns in Table 2 for asymmetric cases are in good agreement between each other. In MD simulations for cases 3 and 4 right after reaching the transition solvophobicity, we see formation of spherical micelles. Upon further increase of solvophobicity in case 4, the rapid coalescence of spherical micelles leads first to formation of worm-like micelles and then macrophase separation (Figure 7). For case 3 the aggregation evolution is similar, except that larger spherical micelles are more stable and result in formation of worm-like micelles with lower degree of shape anisotropy. Such information about the structural evolution of the aggregates is missing in the PRISM theory because the homogeneity basis of the PRISM theory does not allow the application of PRISM theory beyond transition from liquid-like structure.

a The reported errors are from linear fit extrapolation for the data used to obtain the values in first and second columns.

nanoparticles,103 and the deviation between theory and simulations was approximately 15%. In the third column of Table 2 we present the micellization transition calculated from MD simulations through information obtained from the clustering algorithm. The cluster number density plotted in Figure 9b as a function of increasing solvophobicity captures how the clusters/micelles develop and the transition value of solvophobicity that causes the system to achieve stable aggregates. At low εBB values, the ρclust = 0 value describes a fluid disordered state. As εBB increases, the ρclust starts to grow, with the growth starting to occur at different εBB values depending on the copolymer design. After reaching the maximum value, which we call the “cusp”, the ρclust plateaus; the plateau value of ρclust also depends on the copolymer design. We identify the ρclust cusp position with dashed lines in Figure 9b as transition values and list their values in the third column of Table 2. These values are in good agreement with 1/SBB(k*) data from PRISM and MD simulations. The range of εBB where the ρclust decreases from the cusp value to the plateau value corresponds to the solvophobicities that drive some of the clusters to coalesce. While PRISM theory captures these aggregation/micellization transition values well, it is not capable of predicting the higher solvophobicity parameter(s) that can create coalesced or bridged clusters/gels/macrophase-separated structures.

4. CONCLUSION In this paper, we have evaluated the capability of PRISM theory to predict structure and thermodynamics during assembly of linear amphiphilic A−B block copolymers as a function of increasing solvophobicity of the B block, copolymer composition, and sequence via direct comparisons with results from MD simulations. First, we have tested two types of intramolecular correlation functions as input to PRISM theory: intramolecular correlation functions from full-fledged MD simulations (ffMD) and from single chain MD simulations (scMD). Our results showed good agreement between intermolecular distribution functions from either type of input intramolecular correlations, hinting that in these systems J

DOI: 10.1021/acs.macromol.7b01419 Macromolecules XXXX, XXX, XXX−XXX

Macromolecules



the effect from other chains on the intramolecular pair correlation is weak and can be neglected. The PRISM results for real space intermolecular correlations at low solvophobicities showed excellent agreement with MD simulation results. We also observed a good quantitative agreement in reciprocal space structure factors for the same conditions. We have analyzed the behavior of growing microphase structure factor peak at low solvophobicities and connected this behavior with assembled structures observed in MD simulations at high solvophobicities. We have compared PRISM predictions for fluid-microphase/aggregate transition solvophobicity by extrapolating the inverse structure factor value at the microphase peak with transition conditions observed from MD simulations at high solvophobicities and observed good semiquantitative agreement. Using the wide range of copolymer systems that we focused on in this paper, we make the case for PRISM theory as an effective tool for predicting important features of the assembly of amphiphilic block copolymers in solution. (i) PRISM predictions are computationally less expensive than MD simulations, especially because intramolecular correlation information, required as an input to solve PRISM equations, was shown not to need potentially expensive full-fledged MD simulations. (ii) The micro- and macrophase peak position provides information on the type of structure to be expected at high solvophobicities. (iii) The prediction of the solvophobicity required to induce microphase/aggregation transition is within 10% of the MD prediction. (iv) Lastly, PRISM calculations of structure factor may inform about inappropriately small simulation box size to correctly capture large scale behavior. While we only focused on linear copolymers in this paper, for future studies with a wide parameter search, e.g., multiple block sequences, nonlinear architectures, and inclusion of nanoscale additives with varying selectivity to each block of the copolymer, PRISM theory can be a valuable screening tool before conducting more computationally expensive simulations.



REFERENCES

(1) Schweizer, K. S.; Curro, J. G. Integral-equation theory of the structure of polymer melts. Phys. Rev. Lett. 1987, 58 (3), 246. (2) Curro, J. G.; Schweizer, K. S. Theory of polymer melts: an integral equation approach. Macromolecules 1987, 20 (8), 1928−1934. (3) Chandler, D.; Andersen, H. C. Optimized cluster expansions for classical fluids. II. Theory of molecular liquids. J. Chem. Phys. 1972, 57 (5), 1930−1937. (4) Curro, J. G. Intermolecular structure and thermodynamics of vinyl polymer liquids: freely-jointed chains. Macromolecules 1994, 27 (17), 4665−4672. (5) Honnell, K. G.; McCoy, J. D.; Curro, J. G.; Schweizer, K. S.; Narten, A. H.; Habenschuss, A. Local structure of polyethylene melts. J. Chem. Phys. 1991, 94 (6), 4659−4662. (6) Schweizer, K. S.; Curro, J. G. Microscopic theory of the structure, thermodynamics, and apparent χ parameter of polymer blends. Phys. Rev. Lett. 1988, 60 (9), 809. (7) Curro, J. G.; Schweizer, K. S. An integral equation theory of polymer blends: athermal mixtures. Macromolecules 1990, 23 (5), 1402−1411. (8) Rajasekaran, J. J.; Curro, J. G.; Honeycutt, J. Theory for the phase behavior of polyolefin blends: application to the polyethylene/isotactic polypropylene blend. Macromolecules 1995, 28 (20), 6843−6853. (9) Stevenson, C. S.; Curro, J. G.; McCoy, J. D.; Plimpton, S. J. Molecular dynamics simulations of athermal polymer blends: Comparison with integral equation theory. J. Chem. Phys. 1995, 103 (3), 1208−1215. (10) Curro, J. G.; Schweizer, K. S. Integral equation theory for compressible polymer alloys: thermodynamics, scattering, and miscibility of Gaussian chains. Macromolecules 1991, 24 (25), 6736− 6747. (11) Ganesan, V.; Jayaraman, A. Theory and simulation studies of effective interactions, phase behavior and morphology in polymer nanocomposites. Soft Matter 2014, 10 (1), 13−38. (12) David, E. F.; Schweizer, K. S. Integral equation theory of block copolymer liquids. I. General formalism and analytic predictions for symmetric copolymers. J. Chem. Phys. 1994, 100 (10), 7767−7783. (13) David, E. F.; Schweizer, K. S. Influence of conformational asymmetries on local packing and small-angle scattering in athermal diblock copolymer melts. Macromolecules 1995, 28 (11), 3980−3994. (14) David, E. F.; Schweizer, K. S. Solvent-mediated interactions and chain conformations in block copolymer liquids. J. Chem. Soc., Faraday Trans. 1995, 91 (16), 2411−2425. (15) Guenza, M.; Schweizer, K. S. Fluctuations effects in diblock copolymer fluids: Comparison of theories and experiment. J. Chem. Phys. 1997, 106 (17), 7391−7410. (16) Guenza, M.; Schweizer, K. S. Local and microdomain concentration fluctuation effects in block copolymer solutions. Macromolecules 1997, 30 (14), 4205−4219. (17) Chatterjee, A. P.; Schweizer, K. S. Liquid-state theory of semidilute and concentrated polymer solutions. Macromolecules 1998, 31 (7), 2353−2367. (18) Pütz, M.; Curro, J. G.; Grest, G. S. Self-consistent integral equation theory for polyolefins: Comparison to molecular dynamics simulations and x-ray scattering. J. Chem. Phys. 2001, 114 (6), 2847− 2860. (19) Heine, D.; Wu, D. T.; Curro, J. G.; Grest, G. S. Role of intramolecular energy on polyolefin miscibility: Isotactic polypropylene/polyethylene blends. J. Chem. Phys. 2003, 118 (2), 914−924. (20) Curro, J. G.; Schweizer, K. S.; Grest, G. S.; Kremer, K. A comparison between integral equation theory and molecular dynamics simulations of dense, flexible polymer liquids. J. Chem. Phys. 1989, 91 (2), 1357−1364. (21) David, E. F.; Schweizer, K. S. Integral equation theory of block copolymer liquids. II. Numerical results for finite hard-core diameter chains. J. Chem. Phys. 1994, 100 (10), 7784−7795. (22) Gromov, D. G.; de Pablo, J. J. Structure of binary polymer blends: Multiple time step hybrid Monte Carlo simulations and self-

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.macromol.7b01419. Additional details for connectivity criterion parameter determination and additional PRISM and MD results that support our conclusions (PDF)



Article

AUTHOR INFORMATION

Corresponding Author

*E-mail [email protected] (A.J.). ORCID

Daniel J. Beltran-Villegas: 0000-0003-1667-3181 Arthi Jayaraman: 0000-0002-5295-4581 Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS We acknowledge financial support from the National Science Foundation under Grant NSF DMREF-1629156. This work used the HPC Farber supercomputer at the University of Delaware. K

DOI: 10.1021/acs.macromol.7b01419 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules consistent integral-equation theory. J. Chem. Phys. 1995, 103 (18), 8247−8256. (23) Curro, J. G.; Webb, E. B., III; Grest, G. S.; Weinhold, J. D.; Pütz, M.; McCoy, J. D. Comparisons between integral equation theory and molecular dynamics simulations for realistic models of polyethylene liquids. J. Chem. Phys. 1999, 111 (19), 9073−9081. (24) Tillman, P. A.; Rottach, D. R.; McCoy, J. D.; Plimpton, S. J.; Curro, J. G. The structure and thermodynamics of energetically and structurally asymmetric polymer blends. J. Chem. Phys. 1998, 109 (2), 806−814. (25) Tillman, P. A.; Rottach, D. R.; McCoy, J. D.; Plimpton, S. J.; Curro, J. G. The effect of attractions on the structure and thermodynamics of model polymer blends. J. Chem. Phys. 1997, 107 (10), 4024−4032. (26) Yethiraj, A.; Hall, C. K.; Honnell, K. G. Site−site correlations in short chain fluids. J. Chem. Phys. 1990, 93 (6), 4453−4461. (27) Yethiraj, A.; Schweizer, K. S. On the scaling of the critical temperature with the degree of polymerization in symmetric polymer blends. J. Chem. Phys. 1992, 97 (8), 5927−5930. (28) Jayaraman, A.; Nair, N. Integrating PRISM theory and Monte Carlo simulation to study polymer-functionalised particles and polymer nanocomposites. Mol. Simul. 2012, 38 (8−9), 751−761. (29) Martin, T. B.; Jayaraman, A. Using Theory and Simulations To Calculate Effective Interactions in Polymer Nanocomposites with Polymer-Grafted Nanoparticles. Macromolecules 2016, 49 (24), 9684− 9692. (30) Zhao, Q.; Feng, Y.; Zhang, C.; Du, Z.; Tian, M.; Mi, J. Extension of integral equation theory to microphase separation of block copolymers. Mol. Phys. 2015, 113 (8), 880−889. (31) Mai, Y.; Eisenberg, A. Self-assembly of block copolymers. Chem. Soc. Rev. 2012, 41 (18), 5969−5985. (32) Cho, H. K.; Cheong, I. W.; Lee, J. M.; Kim, J. H. Polymeric nanoparticles, micelles and polymersomes from amphiphilic block copolymer. Korean J. Chem. Eng. 2010, 27 (3), 731−740. (33) Letchford, K.; Burt, H. A review of the formation and classification of amphiphilic block copolymer nanoparticulate structures: micelles, nanospheres, nanocapsules and polymersomes. Eur. J. Pharm. Biopharm. 2007, 65 (3), 259−269. (34) Tritschler, U.; Pearce, S.; Gwyther, J.; Whittell, G. R.; Manners, I. 50th Anniversary Perspective: Functional Nanoparticles from the Solution Self-Assembly of Block Copolymers. Macromolecules 2017, 50 (9), 3439−3463. (35) Blanazs, A.; Armes, S. P.; Ryan, A. J. Self-assembled block copolymer aggregates: from micelles to vesicles and their biological applications. Macromol. Rapid Commun. 2009, 30 (4−5), 267−277. (36) Smart, T.; Lomas, H.; Massignani, M.; Flores-Merino, M. V.; Perez, L. R.; Battaglia, G. Block copolymer nanostructures. Nano Today 2008, 3 (3), 38−46. (37) Zhulina, E.; Borisov, O. Theory of block polymer micelles: recent advances and current challenges. Macromolecules 2012, 45 (11), 4429−4440. (38) O’Reilly, R. K.; Hawker, C. J.; Wooley, K. L. Cross-linked block copolymer micelles: functional nanostructures of great potential and versatility. Chem. Soc. Rev. 2006, 35 (11), 1068−1083. (39) Cameron, N. S.; Corbierre, M. K.; Eisenberg, A. 1998 EWR Steacie Award Lecture Asymmetric amphiphilic block copolymers in solution: a morphological wonderland. Can. J. Chem. 1999, 77 (8), 1311−1326. (40) Riess, G. Micellization of block copolymers. Prog. Polym. Sci. 2003, 28 (7), 1107−1170. (41) Raspaud, E.; Lairez, D.; Adam, M.; Carton, J.-P. Triblock copolymers in a selective solvent. 1. Aggregation process in dilute solution. Macromolecules 1994, 27 (11), 2956−2964. (42) Madsen, J.; Armes, S. P. (Meth) acrylic stimulus-responsive block copolymer hydrogels. Soft Matter 2012, 8 (3), 592−605. (43) Kataoka, K.; Harada, A.; Nagasaki, Y. Block copolymer micelles for drug delivery: design, characterization and biological significance. Adv. Drug Delivery Rev. 2001, 47 (1), 113−131.

(44) Allen, C.; Maysinger, D.; Eisenberg, A. Nano-engineering block copolymer aggregates for drug delivery. Colloids Surf., B 1999, 16 (1), 3−27. (45) Jeong, B.; Bae, Y. H.; Lee, D. S.; Kim, S. W. Biodegradable block copolymers as injectable drug-delivery systems. Nature 1997, 388 (6645), 860−862. (46) Chen, S.; Guo, C.; Hu, G.-H.; Liu, H.-Z.; Liang, X.-F.; Wang, J.; Ma, J.-H.; Zheng, L. Dissipative particle dynamics simulation of gold nanoparticles stabilization by PEO−PPO−PEO block copolymer micelles. Colloid Polym. Sci. 2007, 285 (14), 1543−1552. (47) Antonietti, M.; Wenz, E.; Bronstein, L.; Seregina, M. Synthesis and characterization of noble metal colloids in block copolymer micelles. Adv. Mater. 1995, 7 (12), 1000−1005. (48) Schweizer, K.; Curro, J. Integral equation theories of the structure, thermodynamics, and phase transitions of polymer fluids. Adv. Chem. Phys. 1997, 98, 1−142. (49) David, E. F.; Schweizer, K. S. Liquid state theory of thermally driven segregation of conformationally asymmetric diblock copolymer melts. Macromolecules 1997, 30 (17), 5118−5132. (50) Schweizer, K. S.; Curro, J. G. RISM theory of polymer liquids: Analytical results for continuum models of melts and alloys. Chem. Phys. 1990, 149 (1−2), 105−127. (51) Leibler, L. Theory of microphase separation in block copolymers. Macromolecules 1980, 13 (6), 1602−1617. (52) Fredrickson, G. H.; Helfand, E. Fluctuation effects in the theory of microphase separation in block copolymers. J. Chem. Phys. 1987, 87 (1), 697−705. (53) Yethiraj, A.; Schweizer, K. S. Integral equation theory of polymer blends: Numerical investigation of molecular closure approximations. J. Chem. Phys. 1993, 98 (11), 9080−9093. (54) Milchev, A.; Bhattacharya, A.; Binder, K. Formation of block copolymer micelles in solution: A Monte Carlo study of chain length dependence. Macromolecules 2001, 34 (6), 1881−1893. (55) Srinivas, G.; Discher, D. E.; Klein, M. L. Self-assembly and properties of diblock copolymers by coarse-grain molecular dynamics. Nat. Mater. 2004, 3 (9), 638−644. (56) Sheng, Y.-J.; Wang, T.-Y.; Chen, W. M.; Tsao, H.-K. A− B diblock copolymer micelles: Effects of soluble-block length and component compatibility. J. Phys. Chem. B 2007, 111 (37), 10938− 10945. (57) Li, Z.; Dormidontova, E. E. Kinetics of diblock copolymer micellization by dissipative particle dynamics. Macromolecules 2010, 43 (7), 3521−3531. (58) Abbas, S.; Li, Z.; Hassan, H.; Lodge, T. P. Thermoreversible morphology transitions of poly (styrene-b-dimethylsiloxane) diblock copolymer micelles in dilute solution. Macromolecules 2007, 40 (11), 4048−4052. (59) Zhulina, E. B.; Adam, M.; LaRue, I.; Sheiko, S. S.; Rubinstein, M. Diblock copolymer micelles in a dilute solution. Macromolecules 2005, 38 (12), 5330−5351. (60) Mortensen, K. Structural studies of aqueous solutions of PEOPPO-PEO triblock copolymers, their micellar aggregates and mesophases; a small-angle neutron scattering study. J. Phys.: Condens. Matter 1996, 8 (25A), A103. (61) Yang, S.; Yuan, S.; Zhang, X.; Yan, Y. Phase behavior of tri-block copolymers in solution: Mesoscopic simulation study. Colloids Surf., A 2008, 322 (1), 87−96. (62) Liu, C. B.; Gong, C. Y.; Huang, M. J.; Wang, J. W.; Pan, Y. F.; Zhang, Y. D.; Li, G. Z.; Gou, M. L.; Wang, K.; Tu, M. J.; et al. Thermoreversible gel−sol behavior of biodegradable PCL-PEG-PCL triblock copolymer in aqueous solutions. J. Biomed. Mater. Res., Part B 2008, 84 (1), 165−175. (63) Nguyen-Misra, M.; Mattice, W. L. Micellization and gelation of symmetric triblock copolymers with insoluble end blocks. Macromolecules 1995, 28 (5), 1444−1457. (64) Chantawansri, T. L.; Sirk, T. W.; Sliozberg, Y. R. Entangled triblock copolymer gel: Morphological and mechanical properties. J. Chem. Phys. 2013, 138 (2), 024908. L

DOI: 10.1021/acs.macromol.7b01419 Macromolecules XXXX, XXX, XXX−XXX

Article

Macromolecules

(87) Hansen, J.-P.; McDonald, I. R. Theory of Simple Liquids: With Applications to Soft Matter; Academic Press: 2013. (88) Schweizer, K. S.; Yethiraj, A. Polymer reference interaction site model theory: New molecular closures for phase separating fluids and alloys. J. Chem. Phys. 1993, 98 (11), 9053−9079. (89) Singh, C.; Schweizer, K. S.; Yethiraj, A. Fluctuation phenomena in structurally symmetric polymer blends. J. Chem. Phys. 1995, 102 (5), 2187−2208. (90) Hooper, J.; Schweizer, K.; Desai, T.; Koshy, R.; Keblinski, P. Structure, surface excess and effective interactions in polymer nanocomposite melts and concentrated solutions. J. Chem. Phys. 2004, 121 (14), 6986−6997. (91) Hall, L. M.; Anderson, B. J.; Zukoski, C. F.; Schweizer, K. S. Concentration fluctuations, local order, and the collective structure of polymer nanocomposites. Macromolecules 2009, 42 (21), 8435−8442. (92) Martin, T. B.; Jayaraman, A. Polydisperse homopolymer grafts stabilize dispersions of nanoparticles in a chemically identical homopolymer matrix: an integrated theory and simulation study. Soft Matter 2013, 9 (29), 6876−6889. (93) Hindmarsh, A. C.; Brown, P. N.; Grant, K. E.; Lee, S. L.; Serban, R.; Shumaker, D. E.; Woodward, C. S. SUNDIALS: Suite of nonlinear and differential/algebraic equation solvers. ACM Transactions on Mathematical Software (TOMS) 2005, 31 (3), 363−396. (94) Zhou, Y.; Stell, G. Fluids inside a porean integral-equation approach: I. General formalism and hard spheres inside spherical and slit pores. Mol. Phys. 1989, 66 (4), 767−789. (95) Yethiraj, A.; Hall, C. K. Integral equation theory for the adsorption of chain fluids in slitlike pores. J. Chem. Phys. 1991, 95 (5), 3749−3755. (96) Plimpton, S. Fast parallel algorithms for short-range molecular dynamics. J. Comput. Phys. 1995, 117 (1), 1−19. (97) Leach, A. R. Molecular Modelling: Principles and Applications; Pearson Education: 2001. (98) Xu, B.; Li, L.; Yekta, A.; Masoumi, Z.; Kanagalingam, S.; Winnik, M. A.; Zhang, K.; Macdonald, P. M.; Menchen, S. Synthesis, characterization, and rheological behavior of polyethylene glycols end-capped with fluorocarbon hydrophobes. Langmuir 1997, 13 (9), 2447−2456. (99) Newman, M. E.; Ziff, R. M. Fast Monte Carlo algorithm for site or bond percolation. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 2001, 64 (1), 016706. (100) Liu, Z.; Kobayashi, K.; Lodge, T. Dynamic light scattering from asymmetric block copolymers in neutral good solvents: Evidence of association in the disordered state. J. Polym. Sci., Part B: Polym. Phys. 1998, 36 (11), 1831−1837. (101) Bates, F. S. Polymer-polymer phase behavior. Science 1991, 251 (4996), 898−905. (102) Jain, S.; Bates, F. S. On the origins of morphological complexity in block copolymer surfactants. Science 2003, 300 (5618), 460−464. (103) Jayaraman, A.; Schweizer, K. S. Effective interactions, structure, and phase behavior of lightly tethered nanoparticles in polymer melts. Macromolecules 2008, 41 (23), 9430−9438.

(65) Fredrickson, G. H.; Ganesan, V.; Drolet, F. Field-theoretic computer simulation methods for polymers and complex fluids. Macromolecules 2002, 35 (1), 16−39. (66) Huang, C.-I.; Lodge, T. P. Self-consistent calculations of block copolymer solution phase behavior. Macromolecules 1998, 31 (11), 3556−3565. (67) Vavasour, J.; Whitmore, M. Self-consistent field theory of block copolymers with conformational asymmetry. Macromolecules 1993, 26 (25), 7070−7075. (68) Drolet, F.; Fredrickson, G. H. Combinatorial screening of complex block copolymer assembly with self-consistent field theory. Phys. Rev. Lett. 1999, 83 (21), 4317. (69) Zhang, L.; Lin, J.; Lin, S. Self-assembly behavior of amphiphilic block copolymer/nanoparticle mixture in dilute solution studied by self-consistent-field theory/density functional theory. Macromolecules 2007, 40 (15), 5582−5592. (70) Müller, M. Studying amphiphilic self-assembly with soft coarsegrained models. J. Stat. Phys. 2011, 145 (4), 967−1016. (71) Wu, J. Density functional theory for chemical engineering: From capillarity to soft materials. AIChE J. 2006, 52 (3), 1169−1193. (72) Choksi, R.; Ren, X. On the derivation of a density functional theory for microphase separation of diblock copolymers. J. Stat. Phys. 2003, 113 (1), 151−176. (73) Banaszak, M.; Chen, C.; Radosz, M. Copolymer SAFT equation of state. Thermodynamic perturbation theory extended to heterobonded chains. Macromolecules 1996, 29 (20), 6481−6486. (74) Tan, S. P.; Winoto, W.; Radosz, M. Statistical associating fluid theory of homopolymers and block copolymers in compressible solutions: Polystyrene, polybutadiene, polyisoprene, polystyrene-blockpolybutadiene, and polystyrene-block-polyisoprene in propane. J. Phys. Chem. C 2007, 111 (43), 15752−15758. (75) Shukla, K.; Chapman, W. SAFT equation of state for fluid mixtures of hard chain copolymers. Mol. Phys. 1997, 91 (6), 1075− 1082. (76) Gross, J.; Spuhl, O.; Tumakaka, F.; Sadowski, G. Modeling copolymer systems using the perturbed-chain SAFT equation of state. Ind. Eng. Chem. Res. 2003, 42 (6), 1266−1274. (77) Grest, G. S.; Kremer, K. Molecular dynamics simulation for polymers in the presence of a heat bath. Phys. Rev. A: At., Mol., Opt. Phys. 1986, 33 (5), 3628. (78) Ceperley, D.; Kalos, M.; Lebowitz, J. L. Computer simulation of the dynamics of a single polymer chain. Phys. Rev. Lett. 1978, 41 (5), 313. (79) Weeks, J. D.; Chandler, D.; Andersen, H. C. Role of repulsive forces in determining the equilibrium structure of simple liquids. J. Chem. Phys. 1971, 54 (12), 5237−5247. (80) Reddy, G.; Yethiraj, A. Implicit and explicit solvent models for the simulation of dilute polymer solutions. Macromolecules 2006, 39 (24), 8536−8542. (81) Spaeth, J. R.; Kevrekidis, I. G.; Panagiotopoulos, A. Z. A comparison of implicit-and explicit-solvent simulations of self-assembly in block copolymer and solute systems. J. Chem. Phys. 2011, 134 (16), 164902. (82) Schweizer, K.; Curro, J. PRISM theory of the structure, thermodynamics, and phase transitions of polymer liquids and alloys. Atomistic Modeling of Physical Properties 1994, 116, 319−377. (83) Rajasekaran, J. J.; Curro, J. G. Intermolecular packing of freely jointed branched polyalkene melts: a microscopic approach. J. Chem. Soc., Faraday Trans. 1995, 91 (16), 2427−2433. (84) Grayce, C. J.; Yethiraj, A.; Schweizer, K. S. Liquid-state theory of the density dependent conformation of nonpolar linear polymers. J. Chem. Phys. 1994, 100 (9), 6857−6872. (85) Melenkevitz, J.; Schweizer, K.; Curro, J. Self-consistent integral equation theory for the equilibrium properties of polymer solutions. Macromolecules 1993, 26 (23), 6190−6196. (86) Nair, N.; Jayaraman, A. Self-Consistent PRISM Theory− Monte Carlo Simulation Studies of Copolymer Grafted Nanoparticles in a Homopolymer Matrix. Macromolecules 2010, 43 (19), 8251−8263. M

DOI: 10.1021/acs.macromol.7b01419 Macromolecules XXXX, XXX, XXX−XXX