An Efficient and Accurate SCF Algorithm for Block Copolymer Films and Brushes Using Adaptive Discretizations

Self-consistent field (SCF) theory serves as a robust tool for unraveling the intricate behavior exhibited by soft polymeric materials. However, the accuracy and efficiency of SCF calculations are crucially dependent on the numerical methods employed for system discretization and equation-solving. Here, we introduce a simple three dimensional SCF algorithm that uses real-space methods and adaptive discretization, offering improved accuracy and efficiency for simulating polymeric systems at surfaces. Our algorithm’s efficacy is demonstrated through simulations of two distinct polymeric systems, namely, block copolymer (BCP) films and polymer brushes. By enhancing spatial resolution in regions influenced by external forces and employing finer contour discretization at grafting chain ends, we achieve significantly more accurate results at very little additional cost, enabling the study of 3D polymeric systems that were previously computationally challenging. To facilitate the widespread use of the algorithm, we have made our 1D-3D SCF code publicly available.


Introduction
Self-assembled copolymer materials have diverse applications in both industry and daily life due to the wide and complex spectrum of possible morphological patterns into which the copolymer molecules may assemble spontaneously.Self-consistent field (SCF) theory, which employs a simplified chain representation and a mean-field approximation to predict the spatial distribution of polymer segments, has proven to be a reliable and powerful tool for predicting the equilibrium morphology of many polymeric systems [1][2][3].
The SCF calculation iteratively computes space-dependent polymer density and associated potential fields from chain statistics propagators until the self-consistent condition is satisfied [4].The chain propagators satisfy nonlinear modified diffusion equations in the variables space and chain contour ("time"), which have to be numerically solved using spectral or real-space methods.Pure spectral methods based on, e.g., Fourier series-based spectral solutions, have the advantage that they do not require a discretization of the chain contour.Matsen and coworkers have demonstrated their capability in accurately constructing morphology phase diagrams for periodic block-copolymer melts [3,[5][6][7].However, they demand prior knowledge of morphology, assuming symmetry in the considered phase, which limits their applicability for discovering new phases [8][9][10][11].To address these limitations, the pseudo-spectral method was introduced, which does rely on a discretization of the contour and switches between Fourier and real-space representations of the system, utilizing the Fourier representation for the evaluation of gradient terms and the real-space representation for the evaluation of nonlinear self-consistent fields [10,12].In contrast, pure real-space methods discretize the diffusion equation within a simulation box and solve for the solution using finite difference schemes [8,[13][14][15][16].This approach is particularly advantageous for complex polymer systems characterized by symmetry breaking or non-periodic boundary conditions, but can be computational expensive in three dimensions.A comprehensive exploration of the advantages and disadvantages of these methods can be found in Ceniceros and Frederickson's detailed review [11].
Our primary objective in this paper is to enhance the efficiency and accuracy of realspace methods, which heavily rely on the discretization of space/chain contours.The solution of modified diffusion equations in the SCF context typically employs lower-order finite difference; however, unless the discretization is very fine, the accuracy and stability of SCF calculations can suffer, especially when dealing with polymer systems containing sharp interfaces [17,18].For instance, in simulating polymer films, in which the substrate/air interface is often represented by a Dirichlet boundary condition and surface interactions with the polymer are introduced artificially through an external potential field, numerical inaccuracies can significantly affect the calculated free energy of the film [19].Similarly, in the context of polymer brushes [14,20], ensuring proper attachment of the grafting end to the substrate involves fixing an end segment on one space grid using a Dirac delta function as the initial condition for the modified diffusion equations [2,14,21].Achieving convergence to an accurate SCF solution necessitates much finer contour discretization compared to free chains, demanding additional computational effort.
The trade-off between spatial discretization and computational efficiency presents a critical challenge, especially for systems requiring higher-dimensional calculations such as cylinders and spheres in thin copolymer films [19,22,23] or particle-grafted chain polymer brushes with angular or radial-dependent morphologies [21,24,25].To address this challenge, our paper introduces a simple scheme that adaptively increases discretization in the spatial domain where external forces are present and refines the discretization in the contour domain at the grafting point.This approach is similar in spirit to other more sophisticated adaptive methods that have recently been proposed in the literature, such as the use of Oc-Tree data structures [26], polygonal meshes [27], and finite element methods [16,28].By optimizing the spatial resolution according to the system's composition, our approach achieves very high accuracy while keeping computational resources low compared to uniform finite-difference grid methods.In the following sections, we evaluate and demonstrate the effectiveness of our adaptive scheme with two test cases of polymeric systems, namely, block copolymer (BCP) films and polymer brushes.

Background: SCF Equations for Two Test Cases
We use the SCF theory for inhomogeneous systems of Gaussian polymers [29].In the following, we list only the most important equations relevant for our test systems; readers may refer to the literature for the details of the derivations [30].

Test Case 1: Diblock Copolymer Film
We consider an incompressible melt of asymmetric AB diblock copolymer molecules with a degree of polymerization N, which is confined between two flat surfaces.We assume that the majority block A occupies a volume fraction f of each diblock copolymer chain and that both blocks share the same statistical segment length b.In the grand canonical ensemble, the free energy takes the following form [19]: where µ " µ 0 `ln G is the chemical potential, G " ρ c R 3 g is the rescaled dimensionless copolymer density in the bulk, ρ c " n{V is the average molecular number density (with n being the total number of copolymer molecules and V the volume of the film), R g " a Nb 2 {6 is the radius of gyration of the noninteracting copolymer chain, and serves as the spatial length unit throughout the paper, Q is the partition function of a single copolymer chain in the mean field of the other chains, ϕ A prq and ϕ B prq are the local concentrations of the A and B segments at a given point r, and χ is the Flory-Huggins parameter specifying the incompatibility of the two segments.The incompressibility of the BCP melt is ensured by the inverse of the isothermal compressibility parameter κ.The last term on the righthand side of Equation (1) describes the interaction energy with the substrate/interface, with Λ A,B Hprq being the surface field.We assume symmetric boundary wetting conditions; the surface interaction energies with A and B segments are Λ A and Λ B , respectively.The surface field is applied within a depth of ϵ from the two surfaces, where h is the distance between the two surfaces.By finding the extremum of the free energy in Equation ( 1) with respect to ϕ A,B prq, we obtain the fields experienced by the A and B segments: The partition function for a single chain is simply Q " ş dr qpr, sqq : pr, 1 ´sq, with qpr, 1 ´sq and q : pr, sq being the partial partition functions for the first Ns and last Np1 ´sq segments.Note that Q does not depend on the specific choice of s.The propagators qpr, sq and q : pr, 1 ´sq obey the modified diffusion equation

Bqpr, sq
Bs " ∆qpr, sq ´ωprqqpr, sq with initial conditions qpr, 0q " 1 and q : pr, 1q " 1.Here, ∆ represents the Laplacian, r denotes the spatial coordinate in units of R g and 0 ď s ď 1 represents the chain coordinate of the coarse-grained chain segment in the units of the chain contour length L c .The field ωprq " ω A prq for 0 ď s ď f and ωprq " ω B prq for f ă s ď 1.The local concentrations of the A and B segments are simply Starting with an initial guess of the field in Equation ( 5), we first solve the propagator qpr, sq and q : pr, 1 ´sq.Next, the new local concentrations ϕ A prq and ϕ B prq are obtained from Equations ( 6) and ( 7), respectively.Finally, these concentrations are used in Equations ( 3) and (4) to solve for the new fields and the old and new fields are mixed according to a prescription of choice [31].This is repeated until convergence to either a metastable or an equilibrium state.

Test Case 2: Homopolymer Brush
We consider a monodisperse brush solution in which n linear homopolymer chains with a degree of polymerization N are grafted onto a flat substrate at one end.In the canonical ensemble, the free energy of such a system is provided by where ν is the excluded-volume parameter.Similar to the BCP film system, the interactions between the segments and the substrate where the chains are grafted are imposed by an external field Hprq, as in Equation ( 2).The self-consistent equations are determined by the extremum of the free energy, leading to ϕprq " nN Qρ o ż 1 0 ds qpr, sqq : pr, sq, Q " The initial conditions for such a system are provided by q : pr, 0q " δpz ´gp q, where g p is the grafting point and δ is the Dirac delta function.

Adaptive Discretization
We solve the diffusion equation using the semi-implicit Crank-Nicolson discretization scheme combined with the alternating direction implicit method (ADI) [32].For simplicity, we illustrate the method for the one-dimensional diffusion equation.Provided that the discretization points are arranged on a grid, extension of the ADI to three dimensions is straightforward.In one dimension, the diffusion equation reads which can be discretized in s and z as where q m i " qpz i , s m q, with m " t1, Nu being the contour steps in s with contour step size δs m " s m ´sm´1 and i " t1, n z u being the spatial steps in z with spatially varying discretization δz i " z i ´zi´1 .Taking the second-order approximation for the Laplacian ∆ z [33], we have , and c i " We can calculate the propagator q m`1 i at step m `1 from the propagator q m i at step m.The key lies in choosing proper adaptive discretization tailored to the system, which enables gains in both computational time and accuracy.
In the thin film case, where the loss of accuracy mainly arises from the external potential added to mimic interactions with the substrate or air interface, we employ finer discretization only in the z direction where external force is present.We keep a uniform contour discretization where δs{L c " 1{N.Our approach involves adaptive discretization with a total of n z grids, achieved through either a cosine function [16] or a step function.The former has a form of δz i " L z 2 " cosp i´1 n z πq ´cosp i n z πq ı for i " t1, n z u, providing continuous discretization.The latter uses a finer δz s " ϵ α s n z near the surface (z ď ϵ or L z ´ϵ ă z ď L z ) and a coarser δz b " L z ´2ϵ p1´2α s qn z for the rest.The latter provides finer discretization; here, α s represents the fraction of grids allocated to the surface region ϵ.Throughout this work, we refer to these adaptive schemes as "cos" and "step", respectively.
For polymer brushes, additional adaptive contour discretization turns out to be crucial.We utilize δs m " " cosp m´1 2N πq ´cosp m 2N πq ı , m " t1, Nu to discretize the contour.Similar to the approach for thin films, we employ two Dirichlet wall boundaries in the z direction.The separation L z is chosen to be greater than the brush height to ensure that the brush's free end is maintained.The grafting point of the polymer is located at a distance of g p " 0.05 from the substrate.Initially, we employ a uniform grid spacing of δz g " 0.001 for z ď 2g p , utilizing n g " 100 grids to ensure fine discretization at the grafting site.Subsequently, we gradually increase the discretization by 2pn z ´ng q πq ´cosp i 2pn z ´ng q πq ı , i " t1, n z ´ng u.

Block Copolymer Film
We first demonstrate the effectiveness of our adaptive discretization scheme in computing the free energy of BCP films, specifically focusing on lamellar-forming BCP films (with f " 0.5) as an example.Figure 1a-c illustrates the calculated free energy as a function of the film thickness using three different discretization schemes.In the case of uniform discretization, the energy curves diverge for varying discretizations, denoted by different grid numbers.These discrepancies become more pronounced for larger thicknesses due to the coarser discretization.In contrast, our adaptive scheme consistently generates the same free energy curves regardless of the number of discretization grids used.Although a minor discrepancy arises for larger thicknesses in the step-adaptive case, the overall performance remains superior to that of the uniform case.
Figure 1d illustrates the inaccuracies versus discretization, quantified by the shift in energy from the extrapolated energy at spatial discretization δz Ñ 0. To facilitate comparison with the uniform case, the energy shift in the adaptive case is plotted versus the averaged discretization, defined by δz av " L z {n z .In the uniform case, the energy shift can be fitted by a cubic polynomial, which converges to zero with increased discretization.As expected, the accuracy in the energy deteriorates more rapidly with coarser discretization when stronger surface interactions are employed.In contrast, the cos-adaptive scheme consistently produces highly accurate energy values for all tested discretizations, as is evident in the data points overlapping on the zero-error baseline.While surface interaction slightly impacts the shift in energy for larger discretizations, it is inconspicuous compared to the uniform case.
Figure 1e plots the energy shift versus δs{L c for three different discretization schemes and various δz av .The fact that all the curves collapse regardless of the discretization scheme and that δz av indicates that adaptive discretization of the space does not require special treatment of δs.The numerical error caused by δs converges to zero when δs{L c ă 10 ´3 for all three cases.
Our SCF calculations for lamellar BCP films show a significant improvement in accuracy with adaptive discretization, especially when using the cos-adaptive scheme.This allows for coarser discretization, making SCF simulations more efficient and enabling the investigation of multi-layered structures in thicker films.
As a further example, we studied thin films of sphere-forming BCPs.In the bulk, prior SCF studies revealed tiny free energy differences between Hexagonally Close-Packed (HCP) and Face-Centered Cubic (FCC) packings, showing that the HCP phase is the true stable phase [7].Here, we investigate whether this remains true for thin films.Unlike in the case of thin films with lamellar or cylindrical order, the study of sphere packings in thin films requires three-dimensional calculations and large systems, which presents a substantial computational challenge.
Using the cos-adaptive scheme, it is possible to compute the free energy of films containing three layers of spheres (Figure 1a).The close packings of FCC and HCP correspond to ABC and ABA stackings of three layers, as illustrated by the inset drawing in Figure 2a and the SCF-calculated density plot in Figure 1b,c.The free energy curves demonstrate that the HCP phase remains the stable phase in films, i.e., it has the lower free energy.Furthermore, they show that the HCP film has a smaller equilibrium thickness (free energy minimum at h ˚" 9.57 R g ) than the FCC film (h ˚" 9.60 R g ).

Homopolymer Brush
Next, we consider polymer brushes.In this case, inaccuracies in SCF calculations continue to arise due to spatial discretization errors; in addition, the SCF results turn out to suffer from contour discretization errors as well.Specifically, implementing the delta function to graft the chain onto the substrate demands a much smaller δs compared to free polymer chains.Additionally, as the system approaches the strong stretching limit, which is characterized by significant chain interactions and brush thicknesses much larger than the radius of gyration of free chains, R g , SCF calculations typically face convergence issues unless a very fine contour discretization is chosen.
Figure 3a plots the energy shift as a function of the averaged spatial discretization δz av {R g for both uniform and adaptive spatial discretizations and a fixed contour discretization of δs{L c " 0.0001.Similar to the previous example of the thin film, the calculations perform badly in the case of uniform discretization as δz increases.This is particularly evident in scenarios with attractive walls, where more segments are attracted to regions with external potential (the blue curve).Conversely, cos-adaptive δz av consistently yields much more accurate energy calculations, as illustrated by the fact that all filled symbols align with the zero-error baseline.In Figure 3b, we compare the energy shift relative to the average contour discretization δs av for both uniform and cos-adaptive methods; we consider various wall interactions while maintaining a small uniform δz{R g " 0.002 throughout the analysis.For the case of uniform δs, the energy shift converges to zero at around δs{L c " 10 ´5, contrasting with 10 ´3 in the free chain melt case.Moreover, the SCF calculation fails to converge beyond δs{L c Ç 10 ´3 under uniform conditions.However, the adaptive δs scheme is able to accommodate larger discretization values while consistently maintaining significantly higher accuracy for different types of wall interactions, as shown by the empty triangles.
We further extend our calculations to the strong stretching limit.Strong Stretching Theory (SST) predicts a parabolic density profile for the polymer [34,35].The normalized concentration can be written as where ϕ 0 " ?6N 1{2 n bρ o A is the renormalization factor for the concentration and L R g " ´24Nν π 2 ϕ 0

¯1{3
is the absolute thickness of the brush.The SCF-calculated densities are showcased in Figure 3c.In the dilute regime, which is characterized by small L, the SCF captures both the depletion region near the wall and the tail at the top of the brush.As the graft density escalates, the monomer density profile in the brush approaches the prediction of the strong stretching limit.Here, the depletion width narrows, accompanied by a sharp concentration increase from zero.Despite nearing the strong stretching limit, our SCF calculations using the adaptive scheme continue to deliver accurate results.This is evident from the close alignment between the SST theory and the red dashed curve, particularly notable at L " 33.9.

Conclusions
We have proposed a method for incorporating adaptive discretization schemes in real-space SCF calculations to improve both their accuracy and computational efficiency.By implementing finer discretization near surfaces, where strong polymer-substrate interactions occur, the proposed method effectively reduces numerical errors in the calculated free energy.Notably, our study shows that the cos-adaptive scheme consistently achieves superior accuracy, even with spatial discretizations larger than δz av ą 0.1, outperforming the uniform scheme by requiring discretizations that are at least ten times smaller.
To illustrate the potential of our method, we studied the morphologies of sphereforming BCP films, focusing on the question of whether HCP or FCC stacking is more favorable.Our analysis indicates that HCP stacking is the thermodynamically stable state in thin films, at least in the example studied by us (confinement between two attractive surfaces); however, the free energy differences between HCP and FCC stacked films are very small.
When looking at polymer brushes, we again found that adaptive spatial discretization significantly enhances the accuracy of SCF calculations.In addition, adaptive discretization of the contour length parameter turns out to be essential for obtaining accurate results as well as for obtaining converged SCF solutions to begin with when the total number of discretization points is low.Choosing finer values of δs av enables more subtle features to be captured, such as the depletion region near the grafting point, which are easily overlooked when using uniform discretization schemes, especially if δz is comparable to the depletion width.Importantly, our adaptive scheme alleviates the computational burden associated with brush simulations by allowing larger values of the averaged discretization parameter δs av to be chosen, similar to the case of ungrafted chains in the bulk, without sacrificing accuracy.This enhancement in efficiency is particularly noteworthy in view of the fact that the costs of SCF calculations are dominated by the costs of repeatedly solving modified diffusion equations in order to obtain the propagators q and q : .Furthermore, our adaptive scheme remains robust in the strong stretching limit, accommodating scenarios where inter-chain interactions exert substantial influence, resulting in brush heights L " R g .While polymer brush SCF calculations in one dimension remain feasible with manageable computational costs even when using uniform discretization schemes, the integration of adaptive contour discretization and spatial discretization extends the computational capabilities, facilitating SCFT calculations in two or three dimensions for intricate morphologies [21,24,25,36,37].
In summary, our study highlights the crucial role of adaptive discretization schemes in advancing SCF calculations, delivering significant enhancements in accuracy and computational efficiency for various polymeric systems.The versatility of our simple approach is demonstrated by two illustrative examples showcasing its applicability to problems involving interfaces and external potentials.In general, the required accuracy of SCF calculations depends on the specific quantity of interest.In the present paper, we have mainly focused on the free energy, which helps to identify the true equilibrium phase from a set of competing candidate structures.In other applications, accurate prediction of the density profiles of specific chain segments may be more important.
We have introduced two types of adaptive discretization schemes, namely, adaptive spatial discretization and adaptive contour discretization.Adaptive spatial discretization is useful in all situations where density or composition profiles vary strongly within the selected regions in space.On the other hand, adaptive contour discretization is useful in situations where particularly strong variations of the propagator function qpr, sq are expected for well-defined values of s as is the case, for instance, with polymer brushes or copolymers close to junction points that connect different blocks.In general, it would be desirable to couple spatial discretization and contour discretization.This is because the basic equation of any SCF iteration scheme, Equation ( 5), has the form of a modified diffusion equation.In such cases, the contour step should not be chosen independent of the spatial discretization [38,39].In fact, an upper bound for δs should typically scale with δz, as δs " δz 1`ϵ with ϵ " 2 [39].Here, we have chosen a sufficiently small δs, i.e., smaller than the upper bound.In future work, we will explore possibilities for coupling the spatial and contour length discretizations such that the contour discretization is adjusted and becomes finer in those regions of space where the spatial discretization is finer.
It should be noted that the term "adaptive" in the present work refers to situations in which a specific inhomogeneous discretization scheme is chosen at the beginning of an SCF calculation and is not changed thereafter.In fact, we would strongly recommend not proceeding otherwise, as changing the discretization in the middle of an SCF iteration loop is likely to result in convergence problems.However, in algorithms that involve many successive SCF calculations, such as dynamic density functional (DDF) simulations, our scheme could be used to set up dynamically adaptive grids that adjust to the current state of the polymer system.
To foster broader adoption and further refinement, we have made our SCF code publicly available (see below).The code has a modular structure; it provides both adaptive real-space methods and uniform discretization pseudospectral methods, as well as, in the latter case, options to perform dynamic density functional (DDF) simulations following the DDF models used in [40][41][42].It can be applied to arbitrary mixtures of linear multiblock copolymers, and can easily be extended to other polymer architectures as well.

Figure 1 .
Figure 1.(a-c) free energy F{A versus the thickness of the film L z {R g when using the three different discretization schemes described in the text.The data are generated using f " 0.5, µ 0 " 2.55, χN " 20, κN " 25, and Λ A " Λ B " 60.(d) Energy shift per area (∆F{A) as a function of the averaged space discretizations δz av for three different strengths of surface interactions Λ A " Λ B " 0, 60 and 120, where ∆F " Fp∆zq ´Fp0q, with Fp0q being the extrapolated value for N z " 8.The empty circle (˝) and solid triangle (İ) symbols represent uniform discretization and adaptive "cos" discretization, respectively, while the red horizontal line indicates ∆F " 0. (e) Energy shift per area (∆F{A) as a function of contour discretizations δs{L c with Λ A " Λ B " 60.The empty circle (˝), square (˝), and triangle (▽) indicate uniform, cos-adaptive, and step-adaptive discretization, respectively.The solid lines of different colors in (d,e) are polynomial fits of ∆F using f pxq " ax `bx 2 `cx 3 for different surface interaction strengths and δz av .

Figure 2 .
Figure 2. Spherical BCP film consisting of three layers.(a) SCF free energy per area F{A of the film as a function of film thickness L z .The data are generated using f " 0.76, µ 0 " 2.25, χN " 20, κN " 25, Λ A " 60, and Λ B " 5.The optimum thickness, which corresponds to the minimum free energy, is marked by ‹.The inset drawing shows the grid of the FCC (ABC) and HCP (ABA) packing.The green rectangle shows the selected cell for the SCF (b,c): Three-dimensional contour plot of the density of the B-block calculated with the SCF for HCP and FCC packings.

Figure 3 .
Figure3.(a) SCF-calculated energy shift per polymer as a function of the averaged spatial discretization with attractive (Λ " ´100), repulsive (Λ " 100), and neutral (Λ " 0) interactions with the substrate.Here, ϕ 0 " 2 and Nν " 1.(b) SCF-calculated free energy per polymer as a function of the contour discretization.In (a,b), the empty circles correspond to uniform discretization in space (a) and contour (b), while the filled triangles correspond to their counterparts using adaptive discretization.The red horizontal line indicates zero error ∆F " 0. The solid curves in (a,b) are fits to the function f pxq " ax `bx 2 `cx 3 .(c) Normalized segment density ϕpzqL{ϕ 0 R g as a function of z{L (see text for definitions) at Λ " 0. The distance is rescaled by the thickness of the brush estimated by the Strong Stretching Theory (SST).