We develop a multiscale simulation model for diffusion of solutes through porous triblock copolymer membranes. The approach combines two techniques: self-consistent field theory (SCFT) to predict the structure of the self-assembled, solvated membrane and on-lattice kinetic Monte Carlo (kMC) simulations to model diffusion of solutes. Solvation is simulated in SCFT by constraining the glassy membrane matrix while relaxing the brush-like membrane pore coating against the solvent. The kMC simulations capture the resulting solute spatial distribution and concentration-dependent local diffusivity in the polymer-coated pores; we parameterize the latter using particle-based simulations. We apply our approach to simulate solute diffusion through nonequilibrium morphologies of a model triblock copolymer, and we correlate diffusivity with structural descriptors of the morphologies. We also compare the model’s predictions to alternative approaches based on simple lattice random walks and find our multiscale model to be more robust and systematic to parameterize. Our multiscale modeling approach is general and can be readily extended in the future to other chemistries, morphologies, and models for the local solute diffusivity and interactions with the membrane.

Block copolymers self-assemble into microphase-separated ordered structures,1,2 such as hexagonally packed cylinders, that can produce isoporous membranes with higher selectivity and permeability compared to membranes made from homopolymers.3,4 The pore diameter is important for engineering membranes that employ a sieving mechanism for filtration. The membrane pore diameters can be controlled through the composition and molecular weight of the polymer, by incorporating additives, or by using a block copolymer blend.5,6 The pores can also be chemically functionalized, such as with solute-selective ligands, to further improve separation performance.7 As a result, block copolymers are promising materials for fabricating membranes, including mesoscopic membranes used for gas separation8 and ultrafiltration membranes used for water filtration.3,9

There are a number of ways to fabricate block copolymer membranes. One method is to spincoat a block copolymer solution onto a substrate and anneal it, allowing the polymers to relax toward their equilibrium state.10 The final structure can be predicted with knowledge of only the enthalpic interactions (i.e., Flory–Huggins parameters) and the molecular weights of the constituent blocks.11 Pores are then created by etching a sacrificial block, while another block acts as the membrane matrix. Another prominent method combines the highly scalable process of nonsolvent induced phase separation, typically used to form homopolymer membranes, and the self-assembly of block copolymers initiated by the evaporation of the solvent (SNIPS).12–15 Although more scalable, controlling the membrane morphology in a SNIPS process poses a significant challenge because nonequilibrium structures can form, including transient percolation networks or spinodal networks.16 

Diblock copolymers have often been used to make membranes because they are relatively simple to synthesize and have a well-known phase diagram.17 Commonly, one block is chosen to be a glassy material, such as polystyrene, that forms a solid matrix, while the other is a sacrificial material, such as poly(lactic acid), that can be etched to form pores. However, triblock copolymers have recently garnered attention18–22 because the additional architectural and chemical complexity of the polymer offers greater tunability of the membrane. For example, a hydrophilic block, such as poly(ethylene oxide), can be inserted between the glassy matrix block and sacrificial block in order to coat the pores and improve water uptake; additional chemical moieties can also be added to this block to improve selectivity.22 Furthermore, the inclusion of another block in the matrix, e.g., attaching polyisoprene onto poly(styrene-b-4-vinylpyridine) to make poly(isoprene-b-styrene-b-4-vinylpyridine), has been shown to improve mechanical stability.19,21 Triblock copolymers also have an expanded phase space that can lead to advantageous morphologies not seen in diblocks19 and widen the phase window of bicontinuous morphologies with improved toughness.23 These bicontinuous structures also do not require alignment of the domains, unlike cylindrical morphologies. However, as the number of possible morphologies expands, it becomes increasingly necessary to determine which are optimal (or even suitable) for filtration and to understand the relationship between morphology and diffusive transport of solutes through the pores.

Fickian diffusion through cylindrical and lamellar pores is well known to be the one- and two-dimensional equivalent of bulk diffusion, respectively. However, diffusion through real porous media can be complicated by issues such as molecular interactions and confinement.24 Experimentally characterizing diffusion through membrane structures is challenging due to the need to fabricate the structures and measure transport over the required length scales.25 As a result, a variety of simulation studies have been performed to investigate an analogous problem of ion self-diffusion through conducting block copolymers.26–28 Shen et al.26 and Alshammasi and Escobedo27 both simulated the transport of ions through common block copolymer morphologies (lamellar, cylindrical, and gyroid) and found that diffusion through the gyroid network is slower than through oriented lamellae because the increased tortuosity of the network hampers diffusion more than three-dimensional continuity and percolation promotes it. On the other hand, Zhang et al.28 found that once dimensionality was taken into account, the morphology had little influence over the anion diffusivity through block copolymeric ionic liquids and instead depended on the concentration of interfacial anions.

Although these studies provide great insight into the structure–transport relationship for various block copolymer morphologies, they have all relied on either simple random-walk diffusion models that do not capture certain pore-scale effects or more detailed particle-based simulations that cannot easily access the length scales of some membrane morphologies. For example, Howard et al. recently combined self-consistent field theory (SCFT), which was used to generate block copolymer membrane morphologies, with a lattice-based random-walk transport model to study the self-diffusion of a solute.29 To further characterize the important features of each morphology, structural descriptors were calculated and random-forest regression was used to investigate which descriptors had greatest correlation with the self-diffusivity. However, this study neglected the effect of the polymers coating the pores and focused solely on the morphology of the glassy matrix. Pore coatings have been shown to play an important role in water transport through lamellar and cylindrical pores using particle-based dissipative particle dynamics (DPD) simulations;30 however, these simulation techniques are too computationally demanding to apply to more complex membrane morphologies. To this end, it would be highly beneficial to have an efficient high-throughput method for faithfully modeling the diffusion of a solute through various block copolymer structures.

In this study, we develop a multiscale simulation framework, built on our prior work,29 that incorporates the effect of the pore coating on the self-diffusion of solutes through triblock copolymer membranes. Specifically, we focus on nonequilibrium membrane morphologies that are difficult to simulate using more detailed models such as DPD. In order to do this, we present a novel SCFT method for simulating the structural effect of the solvent on the pore coating, then use an on-lattice kinetic Monte Carlo (kMC) model accounting for obstruction from the pore coating to simulate solute diffusion through the pores of the SCFT-generated morphologies. The kMC model uses local pore-level diffusion data from more detailed DPD simulations to model membrane-scale diffusion at much larger length and time scales than are accessible in DPD. Our use of kMC is partially motivated by its recent success in capturing the impact of shale rock features on fluid transport where discrete regimes were used to distinguish between the center of the porous region and the interface.31,32 After determining the self-diffusion coefficient from kMC simulations, we compare the results of our more detailed model with our prior work using simple random-walk simulations based on two definitions of the pore. We show that including the pore coating leads to qualitatively different trends in diffusivity as a function of block fractions, which we correlate with various structural descriptors of the morphology. We find that the features that are important for predicting the solute diffusivity are vastly different when effects of the pore coating are taken into account.

We studied membrane morphologies made from a model ABC triblock copolymer (Fig. 1) having χACN = 35 and χABN = χBCN = 13, where χij is the Flory–Huggins interaction parameter between blocks i and j, and N is the overall degree of polymerization. Tyler et al. computed the equilibrium phase diagram of this idealized nonfrustrated triblock copolymer as a function of the overall block volume fractions fA, fB, and fC using SCFT.33 SCFT refers to a mean-field treatment of a model of interacting polymer molecules in which the degrees of freedom are expressed as auxiliary fields. SCFT is implemented numerically by an algorithm to determine the minimum free energy of the system. Such extremal solutions correspond to saddle points of the Hamiltonian in the space of complex field variables. SCFT yields direct access to free energies and thermodynamic properties, as well as density fields and structure. The majority of phases Tyler et al. discovered were core–shell analogs of the structures found in diblocks, such as lamellae and hexagonally packed cylinders, with one end block (A or C) acting as the core, the middle block (B) acting as the shell, and the other end block (C or A) forming a continuous domain. These structures are of particular interest in membrane fabrication because the core-forming block can be sacrificed, e.g., using etching, to form pores coated by the middle block embedded within a self-supporting continuous matrix. For this study, we designated the A block as a glassy block that forms the membrane matrix and the C block as the sacrificial block to be removed.

FIG. 1.

Multiscale simulation workflow used to study solute diffusion through the solvated pores of the self-assembled ABC triblock copolymer membranes. First, SCFT is used to model self-assembly of the ABC triblock copolymer. The C block is then etched out and replaced with solvent. Another SCFT simulation is run to relax the pore-coating B block and solvent while constraining the glassy A block to maintain the initial matrix morphology. The resulting density distributions are input to the kMC model, along with a model for the local diffusivity computed with DPD simulations, to obtain the diffusivity of solute tracers through the pores.

FIG. 1.

Multiscale simulation workflow used to study solute diffusion through the solvated pores of the self-assembled ABC triblock copolymer membranes. First, SCFT is used to model self-assembly of the ABC triblock copolymer. The C block is then etched out and replaced with solvent. Another SCFT simulation is run to relax the pore-coating B block and solvent while constraining the glassy A block to maintain the initial matrix morphology. The resulting density distributions are input to the kMC model, along with a model for the local diffusivity computed with DPD simulations, to obtain the diffusivity of solute tracers through the pores.

Close modal

In order to study the diffusion of a solute through the porous membrane, we must resolve not only solute-level motion (1 nm) and interactions with pore-level coatings and structures (10 nm), but also consider the effect of the membrane structure as a whole (many times larger than a pore) on diffusion. Since this task is currently intractable for any single simulation method, we developed and applied a multiscale modeling workflow that incorporates various simulation methods that are better-suited to the different scales mentioned (Fig. 1). First, the self-assembly of the triblock copolymer is simulated using SCFT. We then “etch” out the C block and replace it with solvent. In order to obtain the distribution of the solvent and pore-coating B block in the newly formed pores, we run an additional SCFT simulation to reach a locally stable structure while constraining the density profile of the glassy A block so that the matrix morphology does not change. Finally, the density distributions from SCFT and the solute dynamics computed from pore-level DPD simulations are input into a kMC simulation. The solute trajectories resulting from the kMC simulations are used to obtain the self-diffusivity of the solute through the porous network. We describe each of these steps in detail.

Six self-assembled morphologies per triblock composition were selected from the nonequilibrium morphologies generated by Howard et al.29 to facilitate direct comparisons between that study and this one. Details on how these morphologies were generated can be found in their publication.29 We focus here only on morphologies with cubic simulation cells of length 16Rg, where Rg = b(N/6)1/2 is the radius of gyration of an ideal polymer chain and b is the statistical segment length of the triblock copolymer. The simulation cell had periodic boundary conditions in all three dimensions. From the melt morphologies, we removed the C block and replaced it with a solvent (S) that we modeled in SCFT as a “point” polymer, while also choosing N = 100 as the nominal degree of polymerization of the triblock copolymer. After etching the C block, the remaining AB diblock had degree of polymerization NAB = (fA + fB)N. The relative fraction of block i ∈ {A, B} within the polymer was fi/(fA + fB), but the overall block volume fractions remained the same. The Flory–Huggins interaction parameters between the polymer and solvent were chosen to be χAS = 2 and χBS = 0 in order to model a hydrophobic A block and a hydrophilic B block (assuming S is water-like).

To mimic the experimental process of solvating the newly formed membrane pores, SCFT simulations were then performed to relax the B block and solvent to a local free-energy minimum while effectively “freezing” the glassy A block so that the overall membrane morphology did not change. In order to freeze the A block, a novel method was developed and implemented into our in-house SCFT software. A harmonic energy penalty U for the A species was added to the Hamiltonian

βU=κ2ρ0dr(ρA(r)ρA,t(r))2,
(1)

where ρA is the density profile of the A block, ρA,t is the target density profile of the A block (i.e., the output from the original SCFT simulation), and κ determines the strength of the penalty. Moreover, β = 1/(kBT), kB is the Boltzmann constant, T is the temperature, and ρ0 is the total monomer density. For sufficiently large values of κ, the ultimate A density profile becomes quantitatively indistinguishable from that obtained in the original SCFT simulation. To our knowledge, this is the first study that employs such a Hamiltonian to “freeze” a target density profile in order to model the etching and solvation process in SCFT.

The diblock and solvent (AB + S) system was represented using the incompressible multispecies exchange model34 with Gaussian chain statistics, and the exchange-mapped chemical potential fields for the AB + S system were initialized using those obtained from the ABC morphologies. The semi-implicit Siedel scheme35 was used to perform field updates, and the modified diffusion equation was solved using a second-order operator-splitting algorithm36 with contour stepping Δs = 0.01N. A simulation with κ = 600 was run until numerical convergence, and then the result was used to initialize simulations with larger κ. This process of chaining together simulations with increasing κ was continued until the root-mean-square difference between the actual and target density profiles was no greater than 0.01Rg3 at any location in the cell. The final outputs of the SCFT simulations are the local volume fraction profiles ϕi(r) (with i ∈ {A, B, C} before solvation and i ∈ {A, B, S} after solvation) that represent the membrane morphology. These volume fraction profiles are input to the next step of the workflow.

Once the morphologies were generated, we employed a continuum approach to model diffusive motion through the membrane. As a heuristic, we assumed that the solute undergoes Brownian motion, neglecting fast-relaxing inertial motion, and we allowed for a diffusion coefficient that depends on its position. For example, diffusion can be slower in regions of higher polymer concentration where motion is more obstructed. For simplicity, we further assumed that the rate of diffusion was similar in all directions so that solute motion was characterized by a scalar diffusion coefficient D(r) at each position r; this assumption is reasonable for small solutes that are not too close to a surface. To account for influence of the membrane morphology on the solute distribution, we further considered an effective external field ψ(r) that acted on the solutes. We note that D(r) and ψ(r) are not necessarily independent inputs to the model, as both account for effective properties of the solutes that may be correlated, e.g., slower diffusion due to strong attraction. The procedure used to connect the microscale models of D(r) and ψ(r) to the SCFT morphologies is deferred to Sec. II C. Here, we complete the description of the transport model.

Under these assumptions, the probability density ρ(r, t) to find a tracer solute at position r and time t then follows a conservation law (Smoluchowski/Fokker–Planck equation) derived from the Langevin equations of motion37 

ρ(r,t)t=D(r)ρ(r)βψ(r)+D(r)ρ(r).
(2)

The initial condition ρ(r, 0) = δ(rr0) was based on the initial position r0 of the tracer at t = 0, while the boundary conditions were periodic for the morphology, as in the SCFT calculations. At long times, ρ evolves toward a steady state ρ (r) ∼ exp[−βψ(r)], so ψ can be chosen to achieve a targeted equilibrium distribution. Note that Eq. (2) can equivalently be written as

ρ(r,t)t=ρ(r)D(r)ρ(r,t)ρ(r).
(3)

To solve Eq. (2), we simulated the motion of a tracer on a lattice using a kinetic Monte Carlo (kMC) method.38–41 The membrane morphology was discretized into a cubic lattice with edge length , and a tracer was assumed to occupy one lattice site. A tracer starting on site i was allowed to hop along Cartesian axis α to an adjacent site iα± in either the forward (+) or reverse (−) direction with a rate k(iα±|i). The evolution of the probability density ρ(i, t) to find a tracer at site i at time t [the discrete equivalent of ρ(r, t)] is characterized by a master equation for this stochastic process,42,43

ρ(i,t)t=αk(i|iα)ρ(iα,t)+k(i|iα+)ρ(iα+,t)[k(iα|i)+k(iα+|i)]ρ(i,t).
(4)

To choose the move rates, we first imposed detailed balance k(iα±|i)ρ(i)=k(i|iα±)ρ(iα±) using the steady-state distribution ρ. It can be shown that in the limit of small , Eq. (4) then approximates

ρ(i,t)t=ρ(i)αDα(i)2rα2ρ(i,t)ρ(i)
(5)
+vα(i)rαρ(i,t)ρ(i),
(6)

where the derivatives are taken with respect to the α-component of the position coordinate, rα. The diffusivity Dα along direction α is defined in terms of the hopping rates at site i by

Dα(i)=22[k(iα+|i)+k(iα|i)],
(7)

while vα is an effective advection along α defined in terms of the hopping rates at site i by

vα(i)=[k(iα+|i)k(iα|i)].
(8)

Note that because D was assumed to be a scalar, Dα and the sum of the hopping rates at site i must be equal in all directions but vα need not be equal. We compared Eq. (6) to Eq. (3) and chose vα to make the two equivalent. This determined the hopping rates as

k(iα±|i)=D(i)2±12D(i)rαD(i)βψ(i)rα,
(9)

where the partial derivatives are evaluated at i.

Our derivation is general to any external field ψ(r) or scalar diffusivity D(r). The external field enforces the internal membrane morphology (e.g., regions excluded to the solute) that can be determined from measurements or simulations. It can also incorporate interactions with the membrane, such as effective attraction due to chemical functionalization. The diffusivity can be estimated from experiments, empirical diffusion models, or more detailed computer simulations.

In this section, we will describe the specific model solute dynamics that we studied in this work and the methodology used to determine ψ(r) and D(r) for the kMC model described in Sec. II B; however, we emphasize that the framework can be readily extended to incorporate other data sources and other types of solutes. For convenience, we assumed that the solute tracer was chemically similar to the solvent, so its steady-state distribution ρ was directly proportional to ϕS(r) as determined by SCFT; this directly gives βψ(r) = −ln ϕS(r) from the SCFT data. Note that this amounts to a potential that excludes the solute tracer from the walls of the pores and reduces its concentration where the pore-coating is more dense.

To determine D(r), we performed DPD simulations44–46 of diffusion through a single lamellar pore. DPD is a mesoscopic particle-based simulation technique that has been widely used to study block copolymers. Some of us recently used DPD to study pore-level diffusion of water in triblock copolymer membranes,30 revealing that interactions between water and the polymers inside the pore lead to slower local diffusion (i.e., at short times) in regions of higher polymer concentration. This decrease in the local diffusivity leads to a commensurate decrease in the average water diffusion at long times. Here, we used DPD simulations to measure the local tracer diffusivity D(r) that we input to the kMC model, but other molecular modeling approaches could also be applied.

We first constructed a DPD model for the ABC triblock copolymers studied using SCFT. The polymers were represented as linear chains of N = 100 beads with mass m and nominal diameter d connected by springs; each bead was assigned a type (A, B, or C) according to its block. All beads interacted through the standard DPD conservative, random, and dissipative forces.46 The DPD repulsion parameter between beads with the same type was aii = 75 kBT/d, while the DPD repulsion parameter between beads with different types aij was chosen to achieve the desired χijN for the model (see below). The DPD friction parameter for all beads was γi = 4.5 m/τ, where τ=βmd2 is the unit of time. In addition to the standard DPD forces, bonded beads additionally interacted through a harmonic potential ub(r)=k(rr0)2/2 with spring constant k = 100 kBT/d2 and r0 = 1.0 d. All simulations were performed using HOOMD-blue (version 2.6.0) with features extended using azplugins (version 0.8.0).47–49 The integration time step was 0.01 τ, and the bead number density was 3/d3.

In order to map length scales between the DPD model and SCFT calculations, we first prepared a homopolymer melt and measured the radius of gyration Rg and end-to-end distance Re, finding Rg21/2=4.466d and Re21/2=10.92d. Both measurements are consistent with an ideal chain having an effective segment length b = 1.09 d,50 which is slightly larger than the nominal bead diameter. We then carried out the analysis outlined by Groot and Warren to connect aij approximately to χijN for our bead–spring model.46 We performed direct coexistence simulations of two homopolymer oligomers having different bead types, varying the length of the oligomers from 2 to 6 beads and the difference in the repulsion parameter for unlike and like beads Δaij = aijaii in the range 10 ≤ ΔaijN ≤ 30. We initialized equal-sized slabs of each oligomer in an orthorhombic box with square cross section (edge length 10 d) and length 30 d then allowed the mixture to equilibrate for 105τ, during which time the initially separated oligomers partially dissolved in each other. We then sampled configurations every 100 τ during a 105τ production simulation. We computed the average bead density profile with center-of-mass shifting using a bin spacing of 0.5 d,51 extracted the coexistence volume fractions from the bulk region of each slab, and used Flory–Huggins theory to determine χijN.50 As in Groot and Warren’s analysis,46χijN was approximately linear in ΔaijN for all N studied, so we used the best fit of our data χijN = 0.297ΔaijN − 0.124 to choose aij from χijN.

We then created lamellar morphologies of the ABC triblock with fA = 0.5 and 0 ≤ fB ≤ 0.2 in an orthorhombic simulation box with square cross section (edge length 50 d) and using the lamellar spacing computed in the SCFT calculations (about 4.66 Rg or 20.8 d). We first simulated the morphologies for 5 × 104τ. Then, we followed a procedure analogous to the SCFT calculations of freezing the A-block, removing the C-block, and adding solvent. We first “froze” all A beads and any B or C beads that were in the A-rich region of the lamella (defined as being 9d from the center of the C-rich region) by setting their velocities to zero. We then converted all unfrozen C beads to solvent (S) beads, and we removed all bonds between S beads and between S and B beads. The DPD repulsion parameter for the A and S beads was chosen as aAS = 82 kBT/d based on χAS using Groot and Warren’s fit,46 while we used aBS = 75 kBT/d to give χBS ≈ 0. We shifted the center-of-mass velocity of the unfrozen beads to zero and no longer integrated the equations of motion for the frozen beads. We then simulated the unfrozen beads for 5 × 104τ, which allowed the B-block to relax against the solvent as in the SCFT calculations. The volume fraction profiles—computed from configurations sampled every 100 τ during the second half of each simulation using the same procedure as for determining χij—were in excellent agreement between DPD and SCFT for all lamellar morphologies studied, both before [Fig. 2(a)] and after [Fig. 2(b)] solvation. The agreement between the SCFT and DPD volume fraction profiles confirms the interactions of the DPD model are similar to the SCFT model, and so it is reasonable to also use the DPD model to measure solute dynamics.

FIG. 2.

Volume fraction profiles from DPD (points) and SCFT (lines) in lamellar morphologies as a function of position z in the direction normal to the surface (a) before and (b) after solvation when fA = 0.50 and fB = 0.20.

FIG. 2.

Volume fraction profiles from DPD (points) and SCFT (lines) in lamellar morphologies as a function of position z in the direction normal to the surface (a) before and (b) after solvation when fA = 0.50 and fB = 0.20.

Close modal

After preparing the solvated membranes, we measured the local diffusivity D(z) of S beads parallel to the pore surface as a function of position z along the axis normal to the pore surface. Computing the solute-tracer diffusivity from the solvent-bead diffusivity is consistent with our simplifying assumption that the solute is chemically similar to the solvent; an additional bead type could be easily introduced if the solute were chemically different. We computed the parallel mean-squared displacement (MSD) of the S beads Δr2(t)|z0 based on their initial z-position z0 using spatial bins of width 1.0 d. We retained only beads that remained in their initial bin at time t in the average.52 The local diffusivity can be extracted from the time derivative of the MSD dΔr2|z0/dt4D(z0) once the diffusive regime is reached; however, most particles tended to diffuse out of their bins before this point. To address this, we selected a small fraction of the S beads in each bin (75 beads or about 1%) as explicit tracers and tethered their z coordinate to the center of the bin using a harmonic potential with spring constant 64 kBT/d2. The number of tracers and strength of this potential was chosen based on simulations of the bulk solvent so that the number of tracers that remained in their initial bin increased during the relevant measurement window without significantly perturbing the system. We equilibrated the restrained system for 104τ, simulated for another 5 × 104τ and sampled solvent bead configurations every 0.1 τ, then extracted the local diffusivity D(z) from the average value of dΔr2|z0/dt in the time window 500 τ to 1000 τ.

To establish a point of reference for diffusivity in the membrane, we also simulated the bulk diffusivity D0 of the S beads. We equilibrated the solvent in a cubic simulation box with edge length 40 d for 103τ, then sampled configurations every 10 τ during a 5 × 104τ production simulation. We computed the three-dimensional MSD ⟨Δr2(t)⟩ of all S beads, and we extracted the long-time diffusion coefficient from the average of its long-time derivative, d⟨Δr2⟩/dt ∼ 6D0, in the time window 2500 τ to 5000 τ. We will report all values of the diffusivity in the membranes relative to D0.

As in our previous work, the local diffusivity was highest in regions of lower polymer concentration and lowest near the pore surfaces [Fig. 3(a)], which have higher polymer concentration. In prior work, some of us showed that changes in D(z) could be described by treating the B-block in the pore as a Brinkman medium with a mesh size set by the polymer correlation length.30 We found that this picture was unable to fully describe our new simulations, which we suspected was due to the polymer model used here producing rougher pore surfaces that created additional obstructions to the solvent. This surface roughness is captured primarily in the A-block concentration and not the B-block concentration [Fig. 3(b)]. We attempted to modify the Brinkman model to also include the A-block concentration but ultimately found the fit unsatisfactory, possibly due to fundamental differences in the obstructions created by the frozen A-block and dynamic B-block. Accordingly, we posited an empirical model that treated the obstruction from the frozen A-block using a Mackie–Meares-type expression53 and the obstruction from the dynamic B-block using our Brinkman model for an ideal polymer chain,30 

D(ϕA,ϕB)=D01ϕA1+ϕAm1+cϕB+c2ϕB291.
(10)

We fit the parameters m = 1.70 and c = 1.12; m is close to the theoretical exponent of 2 for the standard Mackie–Meares model,53 while c is a fitting parameter accounting for the hydrodynamic radius of the tracer and the scaling prefactor of the polymer correlation length.30,50 This empirical model for D(ϕA, ϕB) was able to fit all our lamellar measurements well [Fig. 3(c)], so we extrapolated it to the more complex membrane morphologies by assuming D(r) = D(ϕA(r), ϕB(r)). We emphasize that this assumption that D depends solely on ϕA and ϕB partially neglects surface curvature effects on D because D was fit for a lamellar morphology. However, it does include curvature to a crude approximation because SCFT does account for curvature when calculating ϕA(r) and ϕB(r). It was not feasible for us to test this approximation, but it is an interesting topic of future investigation.

FIG. 3.

(a) Local diffusivity D(z) in lamellar morphologies with fA = 0.50 and varied fB from DPD simulations. (b) The same data as a function of local A-block volume fraction ϕA for varied fB. (c) The same data as a function of ϕA and local B-block volume fraction ϕB. In (c), the points are colored according to the measured diffusivity and the background shows the fit to Eq. (10). The line in (b) is the fit drawn for ϕB = 0.

FIG. 3.

(a) Local diffusivity D(z) in lamellar morphologies with fA = 0.50 and varied fB from DPD simulations. (b) The same data as a function of local A-block volume fraction ϕA for varied fB. (c) The same data as a function of ϕA and local B-block volume fraction ϕB. In (c), the points are colored according to the measured diffusivity and the background shows the fit to Eq. (10). The line in (b) is the fit drawn for ϕB = 0.

Close modal

We applied the kMC approach of Sec. II B, using the model inputs from Sec. II C to study the diffusion of a solute tracer through the various membrane morphologies prepared in Sec. II A. We discretized the membrane morphology onto the same lattice as was used in the SCFT calculations. (A smaller lattice spacing was tested for selected configurations and found not to significantly impact the results.) We computed ψ using ϕS, and D(r) using ϕA and ϕB at each lattice site. To ensure there was no diffusion into regions of large ϕA, which should be negligible in both experiments and the DPD simulations, we set the external field to βψ = ∞ when ϕS < 0.02 to effectively disallow moves to these sites (make k zero). We then computed the kMC hopping rates according to Eq. (9) using a second-order central finite difference scheme to estimate the required gradients. Care was taken to use an appropriate forward or backward difference when βψ = ∞ at a neighboring site.

To carry out the moves, we employed a rejection-free scheme that randomly selected a move and randomly advanced the time coordinate for each tracer according to the rates at a given lattice site.54–56 For a tracer at lattice site i, we randomly chose to move to an adjacent lattice site iα± with the weight of selecting each site being proportional to k(iα±|i). We unconditionally accepted this move and advanced the time coordinate for each tracer by a random amount that was exponentially distributed with mean [αk(iα±|i)]1, where the sum is over all adjacent sites of i. We found that this rejection-free scheme was computationally advantageous compared to a scheme that advanced the time by a fixed amount but allowed move rejection, i.e., the tracer could remain at site i.56 The rejection scheme was roughly twice as fast as the rejection-free scheme per move because it required half as many random numbers, but the rejection-free scheme was ultimately faster overall when the move rejection rate exceeded 50%. This rejection rate was quite common in the nonequilibrium morphologies, where the average tracer diffusivity was significantly less than the nominal bulk diffusivity D0.

Using this scheme, we simulated an ensemble of 5 × 104 tracers that we randomly initialized onto all lattice sites for which βψ ≠ ∞. We simulated for 104τ to allow the tracers to relax to their equilibrium distribution, then sampled the tracer coordinates every 100 τ during a 105τ simulation. We computed ⟨Δr2⟩ for all tracers and extracted the long-time diffusivity D from d⟨Δr2⟩/dt ∼ 6D in the time window 2 × 104τ to 5 × 104τ. From this point forward, we will distinguish between the long-time diffusivity D and the local diffusivity D(r) based on the inclusion of the position r.

Having developed our multiscale model, we proceeded to analyze the model’s predictions and compare them to alternative approaches. In particular, some of us previously used a simple random-walk (RW) model for solute diffusion in triblock copolymer membranes.29 Unlike the approach described in Sec. II B, the RW model adopted a binary definition of the pores based on the local volume fraction of the A-block matrix ϕA: sites with ϕA < 0.5 (or ϕB + ϕS ≥ 0.5) were defined as the pores, and all other sites were defined as matrix inaccessible to the solute. Within the pores, variations in the distribution of the solute and the local diffusivity were both neglected. This RW model is effectively a special case of the kMC approach assuming βψ = 0 and D/D0 = 1 inside the pores, and βψ = ∞ outside the pores. The RW model is simpler than our new approach because it does not require the additional SCFT calculation to solvate the B-block inside the pore; however, the RW model is potentially less accurate because it neglects pore-level effects on the solute distribution and transport.30 Hence, we compared our diffusion measurements using the new approach that incorporates these effects (Sec. II) to the RW model. We will refer to our new approach as the “kMC model” and the approach of Ref. 29 as the “RW-B+S model;” the latter emphasizes that the pore definition in that model includes both the B block and the solvent.

We first considered the average solute diffusivity D as a function of the polymer block fractions fi (Fig. 4) for the nonequilibrium morphologies of Ref. 29 (Sec. II A). We excluded all morphologies with fA ≤ 0.2 because they tended to produce matrices (using the RW-B + S definition of pore and matrix sites) that were not well connected and so were not mechanically viable for membrane applications. This was not surprising given the limited content of matrix-forming A block and that the equilibrium self-assembled structures for these polymers, such as sphere or disordered phases, also are not connected. We further excluded any morphologies where the lattice sites accessible to the tracer (within each model) were not percolated through the periodic boundaries in at least one dimension, as these morphologies would produce D/D0 = 0 at long times. Figure 4(a) shows the diffusivity predicted by the kMC model, while Fig. 4(b) shows the same for the RW-B + S model. As a verification of our approach for solvating the morphologies in SCFT, we note that we achieved nearly quantitative agreement between Fig. 4(b) and the data of Ref. 29 (largest difference in D/D0 of 0.033). This agreement helps confirm that our method for constraining the matrix-forming A block is effective because simulations of the RW-B + S model using the solvated morphologies should produce identical results to simulations using the original melt morphologies of Ref. 29 if ϕA(r) is successfully frozen.

FIG. 4.

The average diffusivity D of a solute tracer in the solvated morphologies using the (a) new kMC model, (b) RW-B + S model with pore definition ϕB + ϕS ≥ 0.5,29 and (c) RW-S model with pore definition ϕS ≥ 0.5. The block fractions fi are those of the original triblock copolymer. Open circles indicate calculations that were attempted but could not be completed for the model because the lattice sites accessible to the tracer did not percolate for any of the six morphologies considered. The solid lines indicate the equilibrium phase boundaries.33 

FIG. 4.

The average diffusivity D of a solute tracer in the solvated morphologies using the (a) new kMC model, (b) RW-B + S model with pore definition ϕB + ϕS ≥ 0.5,29 and (c) RW-S model with pore definition ϕS ≥ 0.5. The block fractions fi are those of the original triblock copolymer. Open circles indicate calculations that were attempted but could not be completed for the model because the lattice sites accessible to the tracer did not percolate for any of the six morphologies considered. The solid lines indicate the equilibrium phase boundaries.33 

Close modal

In general, D increased in both models as fA decreased due to the increased space available for the tracer to diffuse. However, the diffusivity simulated using the kMC model was smaller than that simulated using the RW-B + S model, particularly for smaller fA, because the local diffusivity D(r) input to the kMC model was typically less than D0 (Fig. 3). To more easily visualize these trends, we projected D as a function of fA along lines of constant fB [Fig. 5(a)], which clearly showed a decrease in D with respect to fA for both models but consistently smaller values of D in the kMC model than in the RW-B+S model.

FIG. 5.

The average diffusivity D of a solute tracer from Fig. 4 projected along lines of constant (a) fB, (b) fA, and (c) fC using the new kMC model, RW-B + S model, and RW-S model. The block fractions fi refer to the original triblock copolymer.

FIG. 5.

The average diffusivity D of a solute tracer from Fig. 4 projected along lines of constant (a) fB, (b) fA, and (c) fC using the new kMC model, RW-B + S model, and RW-S model. The block fractions fi refer to the original triblock copolymer.

Close modal

An even more striking difference between the two models was the dependence of D on fB. Based on prior DPD simulations,30 we expected the diffusivity to decrease with increasing fB because the B block obstructs solute motion in the pores. The kMC model produced diffusivities consistent with this expectation [Fig. 4(a)], but the RW-B+S model did not [Fig. 4(b)]. For example, in Fig. 5(a), the diffusivities obtained using the kMC model were smaller when fB = 0.3 than when fB = 0.1, but the diffusivities obtained using the RW-B+S model were essentially the same. To more clearly highlight this behavior, we projected D as a function of fB along lines of constant fA [Fig. 5(b)]. For the kMC model, there was a clear monotonic decrease in D as a function of fB; however, for the RW-B+S model, there was very little variation in D with fB. In fact, we saw a slight increase in D at the largest values of fB when fA = 0.4. This difference is a consequence of the kMC model incorporating a spatially varying local diffusivity D(r) that depends on the local polymer concentration in the pores rather than the constant diffusivity assumed in the RW-B+S model; we emphasize that we consider the former to be more realistic.

Furthermore, projecting the diffusivity predicted by the kMC model along lines of constant fC revealed that D was roughly constant with respect to changes in fA [Fig. 5(c)]. In order to maintain fixed fC = 1 − fAfB, fB must decrease as fA increases. The local diffusivity in the pores depends on the composition, and the exchange of A for B can lead to a direct competition between the effects of these blocks on the diffusivity. For example, increasing fA tends to decrease D by decreasing the matrix volume that is accessible to the solute; however, the necessary accompanying decrease in fB at constant fC leads to less B in the pores and tends to increase D. In contrast, in the RW-B + S model, the pore size decreases as fA increases at constant fC, but there is no accompanying increase in the diffusivity inside the pores (due to the decrease in fB), resulting in more obstructed motion and smaller D. This compensating behavior may be sensitive to the solute and polymer chemistry (χij) as well as the model for the solute diffusivity.

We also noted that there were some exceptions to this behavior. When fC ≤ 0.2, we no longer observed constant diffusivity along lines of constant fC in the kMC model, primarily when fA was large [Fig. 4(a)]. Increasing fA to large values led to fragmentation of the pores due to the formation of larger A domains, reducing the accessible pathways for diffusion. In this case, exchanging fA for fB at constant fC led to an overall increase in diffusivity. In membrane applications, exchanging the fractions of A and B may also have practical implications on separation performance that are not directly captured by D. For example, if the A block is hydrophobic and the B block is hydrophilic, increasing the B block content may increase water uptake. Moreover, increasing fA may lead to an increase of dead-ends in nonequilibrium morphologies that hamper directed transport across the membrane.57 

After comparing our new approach to our prior work, we asked whether we could combine features of both to construct a RW model that captured similar trends in D as the kMC model but maintained the simplicity of the RW approach. The B-block has two important effects in the kMC model: it excludes the solute from parts of the pore through ψ and it obstructs diffusion in the pore through D(r). The RW model can approximately treat the first effect by redefining the pores. We proposed an alternative model definition where sites with ϕA + ϕB < 0.5 (or ϕS ≥ 0.5) were considered to be the pores, and all other sites were the matrix. This should be considered a crude but convenient approximation because even dense regions of B were typically partially permeable to the solute in all our models. We will refer to this model as the “RW-S” model to emphasize that the pore definition includes essentially only the solvent.

Similar to our results with the kMC and RW-B+S models, the RW-S model captured the monotonic decrease in D as a function of fA [Figs. 4(c) and 5(a)]. Notably, the RW-S model also gave values of D that decreased with increasing fB [Fig. 5(b)], indicating that effects of the B-block can be at least partially approximated as an additional obstruction. However, D was not constant with respect to fA at constant fC [Fig. 5(c)] in the RW-S model. This finding is consistent with our hypothesis that D is nearly constant in the kMC model along lines of constant fC because of a trade-off between the A-block and B-block in the local diffusivity; the RW-S model has constant local diffusivity, so this trade-off cannot be captured. We also note a practical drawback of the RW-S model: because the pore definition is more restrictive than that in the RW-B+S model, we were unable to run simulations for many of the polymers with large fA or fB. Our kMC model, on the other hand, did not suffer from this because an artificial (binary) classification of the lattice sites as pore or matrix was not required. In this respect, we view our new multiscale approach based on the kMC model as being more convenient, systematic, and faithful to physical expectations than either of the RW models.

In Ref. 29, we showed that the diffusivity of the RW-B + S model correlated strongly with two structural descriptors of the pores: namely, their volume v and integrated mean curvature h, normalized by the total membrane volume. These two descriptors, along with the normalized surface area s and integrated Gaussian curvature g, comprise the four Minkowski functionals from integral geometry and image analysis.58,59 Given the qualitative differences we observed between the kMC model and the RW-B+S model, we asked whether the diffusivity of the kMC model correlated with any of these functionals.

The Minkowski functionals can be computed from black-and-white (binary) images using a voxel-counting algorithm.58 The RW-B+S model has a binary pore definition suitable for this algorithm, but the kMC model does not because the solute can access most lattice sites but diffuses at different rates through them. Given the great interest in discovering simple structural descriptors that correlate with diffusivity, we proceeded by adopting an additional structural definition for the pores in the kMC model. Two natural choices were that of either the RW-B+S model or the RW-S model. Given that the RW-S model did not give percolated pores at many state points [Fig. 4(c)], we chose to use the RW-B+S model to define the pores. Hence, the Minkowski functionals we computed were the same for a given morphology for both the kMC and RW-B + S models, but the corresponding diffusivities were different.

We plotted D against the four Minkowski functionals for both models (Fig. 6). In sharp contrast to the RW-B + S model [Figs. 6(e) and 6(f)],29D obtained using the kMC model did not correlate strongly with the volume v or integrated mean curvature h [Figs. 6(a) and 6(b)]. The relationship between D and v was no longer approximately one-to-one: there were many values of v that gave similar values of D. We rationalized this as being due to the B-block coating the pores. Within the RW-B + S model, the pore volume does not change as the length of the pore coating increases and neither does the diffusivity; however, the diffusivity computed using the kMC model does change [Fig. 5(b)]. Similar differences were observed for correlation with h between the two models. Neither the kMC model nor the RW-B+S model showed correlation between D and the surface area s or integrated Gaussian curvature g (not shown).

FIG. 6.

Correlation of the diffusivity D of a solute tracer calculated by the (a)–(d) kMC model and (e)–(f) RW-B + S model with various selected structural descriptors: the pore volume v and integrated mean curvature h (normalized by the total volume of the membrane), the C-block fraction of the triblock copolymer fC, and the average A-block volume fraction experienced by the solute ⟨ϕA⟩. Data are shown for all morphologies that percolated in at least one dimension.

FIG. 6.

Correlation of the diffusivity D of a solute tracer calculated by the (a)–(d) kMC model and (e)–(f) RW-B + S model with various selected structural descriptors: the pore volume v and integrated mean curvature h (normalized by the total volume of the membrane), the C-block fraction of the triblock copolymer fC, and the average A-block volume fraction experienced by the solute ⟨ϕA⟩. Data are shown for all morphologies that percolated in at least one dimension.

Close modal

Motivated by this lack of correlation between D and the Minkowski functionals for the kMC model, we considered additional structural descriptors of the polymer or membrane that might capture variations in D. The simplest descriptors we added were the polymer block fractions fi, which we chose because of the trends observed in Figs. 4 and 5. Only two of these block fractions are independent parameters, so we focus our discussion on fB and fC. The diffusivity was uncorrelated with fB for both models (not shown), which we attribute to the inability of fB to capture major structural changes in the morphology as fA or fC is varied for the triblock copolymer we studied. Interestingly, D correlated well with fC in the kMC model for sufficiently large fC [Fig. 6(c)] although it was essentially uncorrelated with fC in the RW-B+S model [Fig. 6(g)]. We noted, however, that the correlation broke down at smaller values of fC where the pores tended to be less percolated. The block fractions are convenient descriptors based on solely the triblock copolymer architecture, but they do not capture the actual self-assembled membrane morphology or the environment experienced by the solute inside the pores.

Given that the local diffusivity is a function of the local composition [Eq. (10)], we posited two additional descriptors for the solute environment to help improve on the block-fraction descriptors. In particular, we computed the average polymer volume fractions experienced by the solute

ϕi=drϕi(r)ρ(r),
(11)

where ρϕS is the steady state probability distribution of the solute, encoded in the kMC model by ψ. We considered both ⟨ϕA⟩ and ⟨ϕB⟩. For the RW-B + S model, we found that D did not correlate strongly with either ⟨ϕA⟩ [Fig. 6(h)] or ⟨ϕB⟩ (not shown); we somewhat expected this result because the local diffusivity in the RW-B + S model did not depend on the local polymer concentration, so the long-time diffusivity depends primarily on the pore morphology and this dependence is better captured by other descriptors. In the kMC model, D had a strong correlation with ⟨ϕA⟩ [Fig. 6(d)] but not ⟨ϕB⟩ (not shown).

In order to quantify which structural descriptors correlate most strongly with the solute diffusivity for each model, we used the same approach as Ref. 29 and regressed D as a function of the four Minkowski functionals, the three block fractions, and ⟨ϕA⟩ and ⟨ϕB⟩ using a random-forest model.60,61 To determine the importance of each descriptor, we randomized the values of each descriptor and measured the resulting mean decrease in accuracy of the model. The importance of each descriptor for each model (Fig. 7) was determined from the relative magnitude of its mean decrease in accuracy, such that the sum of importances was 1. The results are largely consistent with our visual observations: v and h were the most important descriptors for the RW-B + S model, while fC and ⟨ϕA⟩ were the dominant features for the kMC model. In the kMC model, we found that fC was roughly three times more important than ⟨ϕA⟩, which we interpreted as being a result of fC capturing most of the variation in D except for polymers having sufficiently small fC. It is interesting that the importance of these descriptors changes significantly between the two transport models.

FIG. 7.

Relative importance of each feature in a random forest regression for diffusivity calculated by the kMC model and RW-B + S model.

FIG. 7.

Relative importance of each feature in a random forest regression for diffusivity calculated by the kMC model and RW-B + S model.

Close modal

These results, obtained using our new multiscale model that accounts for pore-level effects on transport, seem to reveal vastly different trends in diffusivity compared to a simpler model that neglected these effects,29 providing new insights into the features that may be important in designing triblock copolymer membranes with nonequilibrium morphologies. We find that the diffusivity of solvent through the pores is mainly impacted by the length of the sacrificial block and the average contact the solvent has with the porous walls. The latter quantity is not readily accessible in experiments, however, it fortunately appears to be most important only when the degree of percolation of the pores is limited. Despite failing to find a direct correlation relationship between diffusivity and structural descriptor based on the Minkowski functionals, we note that fC does ultimately play a role in forming the structure in a way not easily captured by these functionals. Furthermore, we emphasize that our set of Flory–Huggins interaction parameters only roughly approximates one triblock copolymer–solvent system, and different interaction parameters may yield alternative or additional correlations.

We have developed a multiscale model to study solute diffusion through a porous triblock copolymer membrane where the pores have a brush-like coating that interacts with the solute. The method uses SCFT to simulate self-assembly and solvation of the membrane, with the latter implemented using a novel field theory to constrain a target density profile. Then, on-lattice kMC was used to simulate diffusion under the influence of a spatially varying external field and local diffusivity describing the solute interactions with the membrane. We applied our model to simulate solute diffusion through nonequilibrium morphologies of a model ABC triblock copolymer, revealing vastly different trends in diffusivity compared to a simpler geometric model that neglected solute–membrane interaction effects.29 In particular, we found that the diffusivity of solvent through the pores was mainly impacted by the length of the sacrificial polymer block used to create the pores as well as the average contact the solute had with the porous walls. The latter quantity is not readily accessible in experiments; however, it appeared most important when the degree of percolation of the pores was limited.

The model system studied here consists of a nonfrustrated ABC triblock copolymer, a water-like solvent that has strong hydrophobic and hydrophilic interactions with the membrane-matrix A block and pore-coating B block, respectively, and a model solute. The Flory–Huggins interaction parameters χij that we have employed hence correspond to only one possible realization of this system. Our multiscale workflow can be easily applied to study specific polymer, solvent, and solute chemistries with appropriate choice of χij for the SCFT simulations and a model for the local diffusivity, which might be obtained from experimental data or molecular simulations. We also expect that the model could be straightforwardly extended to incorporate other interactions between the solute and membrane, such as electrostatics.

This work was supported as part of the Center for Materials for Water and Energy Systems (M-WET), an Energy Frontier Research Center funded by the U.S. Department of Energy, Office of Science, Basic Energy Sciences under Award No. DE-SC0019272. V.G. and T.M.T. acknowledge financial support from the Welch Foundation (Grant Nos. F-1599 and F-1696). S.K. acknowledges support from the National Science Foundation through the Center for Dynamics and Control of Materials: an NSF Materials Research Science and Engineering Center (NSF MRSEC) under Cooperative Agreement No. DMR-1720595. We acknowledge the Texas Advanced Computing Center (TACC) at The University of Texas at Austin for providing HPC resources. Use was made of computational facilities purchased with funds from the National Science Foundation (Grant No. CNS-1725797) and administered by the Center for Scientific Computing (CSC). The CSC is supported by the California NanoSystems Institute and the Materials Research Science and Engineering Center (MRSEC; NSF Grant No. DMR-1720256) at UC Santa Barbara.

The authors have no conflicts to disclose.

Anthony J. Cooper: Conceptualization (equal); Data curation (lead); Formal analysis (lead); Investigation (equal); Methodology (equal); Visualization (lead); Writing – original draft (equal). Michael P. Howard: Conceptualization (equal); Data curation (supporting); Formal analysis (supporting); Investigation (equal); Methodology (equal); Software (equal); Visualization (supporting); Writing – original draft (equal); Writing – review & editing (equal). Sanket Kadulkar: Methodology (equal); Software (equal); Writing – review & editing (equal). David Zhao: Methodology (supporting); Software (supporting); Writing – review & editing (supporting). Kris T. Delaney: Conceptualization (equal); Methodology (equal); Software (equal); Supervision (equal); Writing – review & editing (equal). Venkat Ganesan: Conceptualization (equal); Funding acquisition (equal); Methodology (equal); Supervision (equal); Writing – review & editing (equal). Thomas M. Truskett: Conceptualization (equal); Funding acquisition (equal); Methodology (equal); Supervision (equal); Writing – review & editing (equal). Glenn H. Fredrickson: Conceptualization (equal); Funding acquisition (equal); Methodology (equal); Supervision (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the authors upon reasonable request.

1.
F. S.
Bates
and
G. H.
Fredrickson
, “
Block copolymer thermodynamics: Theory and experiment
,”
Annu. Rev. Phys. Chem.
41
,
525
557
(
1990
).
2.
F. S.
Bates
, “
Polymer-polymer phase behavior
,”
Science
251
,
898
905
(
1991
).
3.
Y.
Zhang
,
N. E.
Almodovar-Arbelo
,
J. L.
Weidman
,
D. S.
Corti
,
B. W.
Boudouris
, and
W. A.
Phillip
, “
Fit-for-purpose block polymer membranes molecularly engineered for water treatment
,”
npj Clean Water
1
,
2
(
2018
).
4.
V.
Abetz
, “
Isoporous block copolymer membranes
,”
Macromol. Rapid Commun.
36
,
10
22
(
2015
).
5.
Y.
Gu
and
U.
Wiesner
, “
Tailoring pore size of graded mesoporous block copolymer membranes: Moving from ultrafiltration toward nanofiltration
,”
Macromolecules
48
,
6153
6159
(
2015
).
6.
M.
Radjabian
and
V.
Abetz
, “
Tailored pore sizes in integral asymmetric membranes formed by blends of block copolymers
,”
Adv. Mater.
27
,
352
355
(
2015
).
7.
R.
Sujanani
,
M. R.
Landsman
,
S.
Jiao
,
J. D.
Moon
,
M. S.
Shell
,
D. F.
Lawler
,
L. E.
Katz
, and
B. D.
Freeman
, “
Designing solute-tailored selectivity in membranes: Perspectives for water reuse and resource recovery
,”
ACS Macro Lett.
9
,
1709
1717
(
2020
).
8.
M. G.
Buonomenna
,
W.
Yave
, and
G.
Golemme
, “
Some approaches for high performance polymer based membranes for gas separation: Block copolymers, carbon molecular sieves and mixed matrix membranes
,”
RSC Adv.
2
,
10745
10773
(
2012
).
9.
W. A.
Phillip
,
B.
O’Neill
,
M.
Rodwogin
,
M. A.
Hillmyer
, and
E. L.
Cussler
, “
Self-assembled block copolymer thin films as water filtration membranes
,”
ACS Appl. Mater. Interfaces
2
,
847
853
(
2010
).
10.
C.
Sinturel
,
M.
Vayer
,
M.
Morris
, and
M. A.
Hillmyer
, “
Solvent vapor annealing of block polymer thin films
,”
Macromolecules
46
,
5399
5415
(
2013
).
11.
J. D.
Moon
,
B. D.
Freeman
,
C. J.
Hawker
, and
R. A.
Segalman
, “
Can self-assembly address the permeability/selectivity trade-offs in polymer membranes?
,”
Macromolecules
53
,
5649
5654
(
2020
).
12.
K.-V.
Peinemann
,
V.
Abetz
, and
P. F. W.
Simon
, “
Asymmetric superstructure formed in a block copolymer via phase separation
,”
Nat. Mater.
6
,
992
996
(
2007
).
13.
S. P.
Nunes
, “
Block copolymer membranes for aqueous solution applications
,”
Macromolecules
49
,
2905
2916
(
2016
).
14.
Z.
Wang
,
X.
Yao
, and
Y.
Wang
, “
Swelling-induced mesoporous block copolymer membranes with intrinsically active surfaces for size-selective separation
,”
J. Mater. Chem.
22
,
20542
20548
(
2012
).
15.
S.
Woo
,
J.
Kim
,
J.
hyun Lee
, and
J.
Bang
, “
Fabrication of block copolymer membranes via SNIPS process
,”
Korean Chem. Eng. Res.
55
,
214
219
(
2017
).
16.
C.
Stegelmeier
,
V.
Filiz
,
V.
Abetz
,
J.
Perlich
,
A.
Fery
,
P.
Ruckdeschel
,
S.
Rosenfeldt
, and
S.
Förster
, “
Topological paths and transient morphologies during formation of mesoporous block copolymer membranes
,”
Macromolecules
47
,
5566
5577
(
2014
).
17.
E. W.
Cochran
,
C. J.
Garcia-Cervera
, and
G. H.
Fredrickson
, “
Stability of the gyroid phase in diblock copolymers at strong segregation
,”
Macromolecules
39
,
2449
2451
(
2006
).
18.
J.
Rzayev
and
M. A.
Hillmyer
, “
Nanoporous polystyrene containing hydrophilic pores from an ABC triblock copolymer precursor
,”
Macromolecules
38
,
3
5
(
2005
).
19.
W. A.
Phillip
,
R. M.
Dorin
,
J.
Werner
,
E. M. V.
Hoek
,
U.
Wiesner
, and
M.
Elimelech
, “
Tuning structure and properties of graded triblock terpolymer-based mesoporous and hybrid films
,”
Nano Lett.
11
,
2892
2900
(
2011
).
20.
R. M.
Dorin
,
D. S.
Marques
,
H.
Sai
,
U.
Vainio
,
W. A.
Phillip
,
K.-V.
Peinemann
,
S. P.
Nunes
, and
U.
Wiesner
, “
Solution small-angle X-ray scattering as a screening and predictive tool in the fabrication of asymmetric block copolymer membranes
,”
ACS Macro Lett.
1
,
614
617
(
2012
).
21.
Z.
Xu
,
T.
Liu
,
K.
Cao
,
D.
Guo
,
J. M.
Serrano
, and
G.
Liu
, “
Thermally stable and mechanically strong mesoporous films of poly(ether imide)-based triblock copolymers
,”
ACS Appl. Polym. Mater.
2
,
1398
1405
(
2020
).
22.
Y.
Zhang
,
R. A.
Mulvenna
,
S.
Qu
,
B. W.
Boudouris
, and
W. A.
Phillip
, “
Block polymer membranes functionalized with nanoconfined polyelectrolyte brushes achieve sub-nanometer selectivity
,”
ACS Macro Lett.
6
,
726
732
(
2017
).
23.
B. J.
Dair
,
C. C.
Honeker
,
D. B.
Alward
,
A.
Avgeropoulos
,
N.
Hadjichristidis
,
L. J.
Fetters
,
M.
Capel
, and
E. L.
Thomas
, “
Mechanical properties and deformation behavior of the double gyroid phase in unoriented thermoplastic elastomers
,”
Macromolecules
32
,
8145
8152
(
1999
).
24.
D. M.
Tartakovsky
and
M.
Dentz
, “
Diffusion in porous media: Phenomena and mechanisms
,”
Transp. Porous Media
130
,
105
127
(
2019
).
25.
D.
Zalami
,
O.
Grimm
,
F. H.
Schacher
,
U.
Gerken
, and
J.
Köhler
, “
Non-invasive study of the three-dimensional structure of nanoporous triblock terpolymer membranes
,”
Soft Matter
14
,
9750
9754
(
2018
).
26.
K.-H.
Shen
,
J. R.
Brown
, and
L. M.
Hall
, “
Diffusion in lamellae, cylinders, and double gyroid block copolymer nanostructures
,”
ACS Macro Lett.
7
,
1092
1098
(
2018
).
27.
M. S.
Alshammasi
and
F. A.
Escobedo
, “
Correlation between ionic mobility and microstructure in block copolymers. A coarse-grained modeling study
,”
Macromolecules
51
,
9213
9221
(
2018
).
28.
Z.
Zhang
,
J.
Krajniak
, and
V.
Ganesan
, “
A multiscale simulation study of influence of morphology on ion transport in block copolymeric ionic liquids
,”
Macromolecules
54
,
4997
5010
(
2021
).
29.
M. P.
Howard
,
J.
Lequieu
,
K. T.
Delaney
,
V.
Ganesan
,
G. H.
Fredrickson
, and
T. M.
Truskett
, “
Connecting solute diffusion to morphology in triblock copolymer membranes
,”
Macromolecules
53
,
2336
2343
(
2020
).
30.
D.
Aryal
,
M. P.
Howard
,
R.
Samanta
,
S.
Antoine
,
R.
Segalman
,
T. M.
Truskett
, and
V.
Ganesan
, “
Influence of pore morphology on the diffusion of water in triblock copolymer membranes
,”
J. Chem. Phys.
152
,
014904
(
2020
).
31.
M.
Apostolopoulou
,
M. S.
Santos
,
M.
Hamza
,
T.
Bui
,
I. G.
Economou
,
M.
Stamatakis
, and
A.
Striolo
, “
Quantifying pore width effects on diffusivity via a novel 3D stochastic approach with input from atomistic molecular dynamics simulations
,”
J. Chem. Theory Comput.
15
,
6907
6922
(
2019
).
32.
M.
Apostolopoulou
,
M.
Stamatakis
,
A.
Striolo
,
R.
Dusterhoft
,
R.
Hull
, and
R.
Day
, “
A novel modeling approach to stochastically evaluate the impact of pore network geometry, chemistry and topology on fluid transport
,”
Transp. Porous Media
136
,
495
520
(
2021
).
33.
C. A.
Tyler
,
J.
Qin
,
F. S.
Bates
, and
D. C.
Morse
, “
SCFT study of nonfrustrated abc triblock copolymer melts
,”
Macromolecules
40
,
4654
4668
(
2007
).
34.
D.
Düchs
,
K. T.
Delaney
, and
G. H.
Fredrickson
, “
A multi-species exchange model for fully fluctuating polymer field theory simulations
,”
J. Chem. Phys.
141
,
174103
(
2014
).
35.
H. D.
Ceniceros
and
G. H.
Fredrickson
, “
Numerical solution of polymer self-consistent field theory
,”
Multiscale Model. Simul.
2
,
452
474
(
2004
).
36.
K. Ø.
Rasmussen
and
G.
Kalosakas
, “
Improved numerical algorithm for exploring block copolymer mesophases
,”
J. Polym. Sci., Part B: Polym. Phys.
40
,
1777
1783
(
2002
).
37.
M. P.
Howard
,
A.
Gautam
,
A. Z.
Panagiotopoulos
, and
A.
Nikoubashman
, “
Axial dispersion of Brownian colloids in microfluidic channels
,”
Phys. Rev. Fluids
1
,
044203
(
2016
).
38.
A. La
Magna
,
S.
Coffa
, and
L.
Colombo
, “
A lattice kinetic Monte Carlo code for the description of vacancy diffusion and self-organization in Si
,”
Nucl. Instrum. Methods Phys. Res., Sect. B
148
,
262
267
(
1999
).
39.
U.
Von Toussaint
,
T.
Schwarz-Selinger
, and
K.
Schmid
, “
First-passage kinetic Monte Carlo on lattices: Hydrogen transport in lattices with traps
,”
J. Nucl. Mater.
463
,
1075
1079
(
2015
).
40.
S.
Kadulkar
,
D.
Banerjee
,
F.
Khabaz
,
R. T.
Bonnecaze
,
T. M.
Truskett
, and
V.
Ganesan
, “
Influence of morphology of colloidal nanoparticle gels on ion transport and rheology
,”
J. Chem. Phys.
150
,
214903
(
2019
).
41.
S.
Kadulkar
,
D. J.
Milliron
,
T. M.
Truskett
, and
V.
Ganesan
, “
Transport mechanisms underlying ionic conductivity in nanoparticle-based single-ion electrolytes
,”
J. Phys. Chem. Lett.
11
,
6970
6975
(
2020
).
42.
V.
Ruiz Barlett
,
M.
Hoyuelos
, and
H. O.
Mártin
, “
Comparison between fixed and Gaussian step length in Monte Carlo simulations for diffusion processes
,”
J. Comput. Phys.
230
,
3719
3726
(
2011
).
43.
V.
Ruiz Barlett
,
M.
Hoyuelos
, and
H. O.
Mártin
, “
Monte Carlo simulation with fixed step length for diffusion processes in nonhomogeneous media
,”
J. Comput. Phys.
239
,
51
56
(
2013
).
44.
P. J.
Hoogerbrugge
and
J. M. V. A.
Koelman
, “
Simulating microscopic hydrodynamic phenomena with dissipative particle dynamics
,”
Europhys. Lett.
19
,
155
160
(
1992
).
45.
P.
Español
and
P.
Warren
, “
Statistical mechanics of dissipative particle dynamics
,”
Europhys. Lett.
30
,
191
196
(
1995
).
46.
R. D.
Groot
and
P. B.
Warren
, “
Dissipative particle dynamics: Bridging the gap between atomistic and mesoscopic simulation
,”
J. Chem. Phys.
107
,
4423
4435
(
1997
).
47.
J. A.
Anderson
,
J.
Glaser
, and
S. C.
Glotzer
, “
HOOMD-blue: A python package for high-performance molecular dynamics and hard particle Monte Carlo simulations
,”
Comput. Mater. Sci.
173
,
109363
(
2020
).
48.
C. L.
Phillips
,
J. A.
Anderson
, and
S. C.
Glotzer
, “
Pseudo-random number generation for brownian dynamics and dissipative particle dynamics simulations on GPU devices
,”
J. Comput. Phys.
230
,
7191
7201
(
2011
).
49.
Dataset: (2022). “A HOOMD-blue component for soft matter simulations,” Github. https://github.com/mphowardlab/azplugins.
50.
M.
Rubinstein
and
R. H.
Colby
,
Polymer Physics
(
Oxford University Press
,
2003
).
51.
K. S.
Silmore
,
M. P.
Howard
, and
A. Z.
Panagiotopoulos
, “
Vapour–liquid phase equilibrium and surface tension of fully flexible Lennard-Jones chains
,”
Mol. Phys.
115
,
320
327
(
2017
).
52.
P.
Liu
,
E.
Harder
, and
B. J.
Berne
, “
On the calculation of diffusion coefficients in confined fluids and interfaces with an application to the liquid–vapor interface of water
,”
J. Phys. Chem. B
108
,
6595
6602
(
2004
).
53.
J. S.
Mackie
and
P.
Meares
, “
The diffusion of electrolytes in a cation-exchange resin membrane. I. Theoretical
,”
Proc. R. Soc. London, Ser. A
232
,
498
509
(
1955
).
54.
D. T.
Gillespie
, “
A general method for numerically simulating the stochastic time evolution of coupled chemical reactions
,”
J. Comput. Phys.
22
,
403
434
(
1976
).
55.
D. T.
Gillespie
, “
Exact stochastic simulation of coupled chemical reactions
,”
J. Phys. Chem.
81
,
2340
2361
(
1977
).
56.
A.
Chatterjee
and
D. G.
Vlachos
, “
An overview of spatial microscopic and accelerated kinetic Monte Carlo methods
,”
J. Comput.-Aided Mater. Des.
14
,
253
308
(
2007
).
57.
L. Y.
Schneider
and
M.
Müller
, “
Engineering scale simulation of nonequilibrium network phases for battery electrolytes
,”
Macromolecules
52
,
2050
2062
(
2019
).
58.
K.
Michielsen
and
H.
De Raedt
, “
Integral-geometry morphological image analysis
,”
Phys. Rep.
347
,
461
538
(
2001
).
59.
R. T.
Armstrong
,
J. E.
McClure
,
V.
Robins
,
Z.
Liu
,
C. H.
Arns
,
S.
Schlüter
, and
S.
Berg
, “
Porous media characterization using Minkowski functionals: Theories, applications and future directions
,”
Transp. Porous Media
130
,
305
335
(
2019
).
60.
L.
Breiman
, “
Random forests
,”
Mach. Learn.
45
,
5
32
(
2001
).
61.
F.
Pedregosa
,
G.
Varoquaux
,
A.
Gramfort
,
V.
Michel
,
B.
Thirion
,
O.
Grisel
,
M.
Blondel
,
P.
Prettenhofer
,
R.
Weiss
,
V.
Dubourg
,
J.
Vanderplas
,
A.
Passos
,
D.
Cournapeau
,
M.
Brucher
,
M.
Perrot
, and
E.
Duchesnay
, “
Scikit-learn: Machine learning in python
,”
J. Mach. Learn. Res.
12
,
2825
2830
(
2011
).