# RuSseL: A Self-Consistent Field Theory Code for Inhomogeneous Polymer Interphases

^{*}

^{†}

Previous Article in Journal

School of Chemical Engineering, National Technical University of Athens (NTUA), GR-15780 Athens, Greece

Authors to whom correspondence should be addressed.

These authors contributed equally to this work.

Academic Editor: Brendan Howlin

Received: 9 April 2021 / Revised: 4 May 2021 / Accepted: 7 May 2021 / Published: 10 May 2021

(This article belongs to the Section Computational Chemistry)

In this article, we publish the one-dimensional version of our in-house code, RuSseL, which has been developed to address polymeric interfaces through Self-Consistent Field calculations. RuSseL can be used for a wide variety of systems in planar and spherical geometries, such as free films, cavities, adsorbed polymer films, polymer-grafted surfaces, and nanoparticles in melt and vacuum phases. The code includes a wide variety of functional potentials for the description of solid–polymer interactions, allowing the user to tune the density profiles and the degree of wetting by the polymer melt. Based on the solution of the Edwards diffusion equation, the equilibrium structural properties and thermodynamics of polymer melts in contact with solid or gas surfaces can be described. We have extended the formulation of Schmid to investigate systems comprising polymer chains, which are chemically grafted on the solid surfaces. We present important details concerning the iterative scheme required to equilibrate the self-consistent field and provide a thorough description of the code. This article will serve as a technical reference for our works addressing one-dimensional polymer interphases with Self-Consistent Field theory. It has been prepared as a guide to anyone who wishes to reproduce our calculations. To this end, we discuss the current possibilities of the code, its performance, and some thoughts for future extensions.

Mesoscopic simulations have been instrumental in unveiling the phenomena and mechanisms that govern the thermodynamics and dynamics of polymer interfaces [1,2]. Over the years, a wide class of such mesoscopic approaches have been developed for addressing such systems, including Lattice Boltzmann [3], particle (e.g., Dissipative Particle Dynamics [4]), hybrid particle-field [5,6,7,8], and purely field theoretic methods, such as density functional theory [9] and Self-Consistent Field Theory (SCFT) [10,11,12,13]. Regarding the latter, mesoscopic calculations of polymer systems based on SCFT have been established as a powerful theoretical tool for the description of their thermodynamic and structural properties [12,14,15]. When applied in dense enough systems, SCFT describes quite accurately the entropic and enthalpic phenomena taking place in complex polymer systems. This is the reason why it has been widely used to address the equilibrium behavior of polymer melts and block copolymers, which are important for industrial applications such as microelectronics and lithography [16,17,18,19]. Furthermore, a considerable number of computational studies based on SCFT have addressed heterogeneous systems involving polymer–solid or polymer–gas interfaces [10,11,12,13].

It is quite usual to formulate SCFT in Fourier space by invoking the so-called spectral and pseudo-spectral methods [20,21]. Applying a Fast-Fourier Transform (FFT) numerical scheme, one can obtain considerable performance in the solution of the partial differential equation (PDE) describing the propagation of chain conformations [21]. Nevertheless, these spectral methodologies require specific boundary conditions (BCs) and geometric symmetries to achieve optimal efficiency [19,22]. Real-space methods, on the other hand, such as Finite Differences (FD) or the Finite Element Method (FEM) [23,24], may be more complex to implement, but offer higher flexibility in the types of physical systems that can be addressed. This is why real-space methods were adopted in this work. The three-dimensional version of RuSseL itself is based on the FEM [25]. Daoulas and Müller have also used a compressible Self-Consistent Field formulation implemented with a three-dimensional finite differences scheme in order to investigate the thermodynamic stability of amphiphilic membranes [26]. Ouaknin et al. have implemented a hybrid finite difference/volume scheme to address block copolymer systems in confined domains, imposing Robin or effective Neumann boundary conditions [15].

The present article publishes the in-house developed code, RuSseL, for field theoretic calculations on inhomogeneous polymer interphases and discusses some of the new functionalities, for example, the determination of the formation of polymer bridges inside the system or the use of tabulated potentials for the description of solid–polymer interactions. The article provides a detailed description of the code, provides several examples and applications in a broad variety of polymeric interfaces, and is supported by a detailed manual. The code is named after the mathematician Bertrand Russell and has already been employed by the authors to investigate a wide range of systems such as polymer–vacuum [11] and polymer–solid interfaces in planar and spherical geometries, polymer-grafted surfaces and nanoparticles (NPs) [27], and systems with opposing surfaces either bare or grafted, embedded in melt or vacuum phases [28]. It should be stated that the smearing approximation adopted for grafting points might not be suitable in the limit of low grafting densities, albeit the qualitative behavior is not expected to vary. Serious efforts were conducted by the authors to investigate these effects with the three dimensional version of RuSseL, which is also going to be published in the future [25].

It contrast to other SCFT codes in the literature [29,30], which, in most cases, evaluate the polymer thermodynamics under the incompressibility condition, RuSseL has been designed and optimized to model any kind of compressible thermoplastic melt with realistic isothermal compressibilities, given that its volumetric properties are provided by the user and/or calculated via an appropriate equation of state.

Even though the compressibility condition makes the convergence of SCFT significantly harder, it offers the ability to address systems comprising gaseous phases, such as those of isolated nanoparticles and grafted surfaces in the absence of a melt phase. These systems are essential to investigating and understanding the structure of particle solids [31] and dry brushes, and to quantifying the solvation free energy of embedded grafted nanoparticles. At the moment, the code supports the Helfand (HFD) and Sanchez–Lacombe (SL) equations of state to describe nonbonded interactions, but it can be easily extended to adopt any equation of state given that its free energy density is provided.

In single surfaces, the code can predict the interfacial free energies, such as surface tension [11] and adhesion tension [28], as functions of the curvature. These properties, when evaluated in conjunction with square gradient theory (SGT), have been shown [11] to be in excellent agreement with experimental measurements by using a single value for the influence parameter. In addition, having calculated the surface or adhesion tension, one can estimate seamlessly the macroscopic wetting functions (i.e., work of cohesion, adhesion, spreading, and immersion), as well as the contact angles [28].

Recently, the code has been employed to quantify the potential of mean force for opposing surfaces, either bare or grafted, embedded in melt or vacuum phases, as a function of the length of the matrix and grafted chains, the grafting densities of the latter, the asymmetry in grafting densities and chain lengths between the two surfaces, and the solvent and wetting conditions [28].

Besides the thermodynamic calculations of interfacial free energies, the code exports several measures concerning the chain structure such as the following:

- density profiles of each chain type
- density profiles of specific chain segments; that is, chain ends, or any other segment specified by the user

The code has been equipped with several kinds of wall potentials, such as the Hamaker potential [27,36], simpler potentials, such as the square well [10,32] and the ramp [28], as well as combinations of the above potentials (hybrid) [28]. Moreover, it offers the possibility to use custom—user-designed—potentials, or tabulated potentials.

The article is structured as follows: Section 2 discusses representative models that can be described using the code, Section 3 presents briefly the underlying mathematical formulation, Section 4 illustrates the overall structure and flow of the code, Section 5 provides information regarding the computes, and Section 6 concludes the manuscript.

At the moment, the code can support three kinds of chains (c) at the same time: matrix chains (c = m), and grafted chains at the low/high boundary (c = ${\mathrm{g}}^{\mp}$). The lengths of the matrix (N_{m}) and grafted (${N}_{{\mathrm{g}}^{\mp}}$) chains, as well as the corresponding grafting density at each surface (${\sigma}_{{\mathrm{g}}^{\mp}}$), can be tuned independently of each other. The type of the interface (i.e., gas vs. solid) is set by choosing the appropriate boundary conditions for each chain kind at the interfaces. Dirichlet boundary conditions, q_{c}(**r**, N) = 0, are applicable for polymer–vacuum and polymer–solid interfaces, where $r\in \partial \Omega $. Neumann boundary conditions, ${\nabla}_{r}{q}_{c}\left(r,N\right)=0$, are applicable for boundaries corresponding to bulk phases (of melt or gas). In addition, the user can choose between planar and spherical geometries, corresponding to planar interfaces or films and nanoparticles or cavities, respectively. The flexibility of tuning these parameters individually enhances the accessible model space considerably.

To facilitate the discussion, the different kinds of boundaries and domains will be denoted as V: vacuum, S: solid surface, G: grafted (solid) surface, and M: melt or matrix; whereas combinations of these will be used to refer to each system. Figure 1 depicts illustrations and density profiles of representative planar geometries with single and opposing surfaces. Similarly, Figure 2 depicts some representative systems for spherical geometries (cavities and nanoparticles).

Setting vacuum interfaces at one or both sides of the domain allows for the modelling of vacuum–matrix interphases (Figure 2, VM) or free films (Figure 2, VMV), respectively [10,11,33,37]. In planar geometries, besides the density profiles and other relevant structural properties, one can estimate the interfacial free energy, that is, the surface tension, in this case. In spherical geometries (Figure 2, VM), the system corresponds to a cavity of radius R embedded in a melt, from which one can estimate the curvature-dependent structural properties of the polymer and the corresponding surface tension; hence, one can obtain a measure regarding the stability of the cavity.

On the other hand, introducing the polymer–solid interactions allows for the modelling of solid–polymer interphases, for example, see SM in Figure 1 and Figure 2. The wettability of these interphases can be tuned by adjusting the kind and the strength of the polymer–solid interactions; we elaborate on this subject in the next section. Enabling the solid interactions on both sides (Figure 1, SMS), allows for the estimation of the potential of mean force (or the disjoining pressure) and thus, obtaining a measure of the stability of such systems [28,34,38].

Introducing grafted chains on one surface (Figure 1 and Figure 2, GM) allows for the modeling of dilute grafted nanoparticles as a function of ${\sigma}_{{\mathrm{g}}^{\mp}}$, ${N}_{{\mathrm{g}}^{\mp}}$, N_{m}, and curvature, and extracts several key structural (density profiles, brush thickness) and thermodynamic quantities [27]. By grafting both surfaces (Figure 1, GMG), one can estimate the potential of mean force as a function of the surface–surface distance, the chain lengths and grafting densities of the grafted chains, and the chain lengths of matrix chains.

Removing the melt phase entirely, we can model grafted surfaces or NPs in dry solvent conditions (Figure 1 and Figure 2, GV). Having calculated the interfacial energies of GM and GV systems, one can estimate the solvation free energy. GVG systems, on the other hand, can provide information about the structure of the so-called “particle solids” (in cases where the particles are extremely coarse) [31,39].

The user is, of course, not limited to choosing these combinations. For example, one might want to examine SMV [37] or even GMV systems.

In addition to the kind of boundary conditions, the code offers the possibility of selecting from a range of different wall potentials. In this way, the wetting degree of the interface can be altered according to the properties of the specific system that the user wishes to reproduce. The current section presents evaluations with the ramp, the square well, and tabulated potentials for demonstration purposes, but details concerning other functional forms, such as the Hamaker potential are also provided in the Supporting Information.

Figure 3 illustrates the density profiles of polyethylene (PE) melt near solid surfaces as a function of the contact angle using the HFD, SL, and SL-SGT EoS. The polymer–solid interactions have been described with the square well (v_{square_well}, left) or the ramp (v_{ramp}, right) potential. The functional forms of the square well and ramp potentials are the following:
where v_{square_well/}v_{ramp} and σ_{square_well/}σ_{ramp} correspond to the depth and the width of the square well/ramp potential. The range of the potentials was set to 6.5 Å, and the position of the hard sphere wall at 2 Å (thus, the minimal distance of a segment from the hard sphere wall is equal to 4.5 Å) [10,32]. For each one of these cases, the depth of the potential was optimized with the Secant optimization method in order to match the target contact angles, $\theta =\mathrm{acos}\left(-{\gamma}^{\mathrm{SM}}/{\gamma}^{\mathrm{VM}}\right)$, with γ^{SM} and γ^{VM} being the interfacial free energy of the SM and VM system, respectively. Note that, γ^{VM} ≡ σ^{VM} is the surface tension, and γ^{SM} ≡ –(σ^{SV} − σ^{SM}) [40] corresponds to minus the adhesion tension; thus, the contact angle can be estimated equivalently as $\theta =\mathrm{acos}\left[\left({\sigma}^{\mathrm{SV}}-{\sigma}^{\mathrm{SM}}\right)/{\sigma}^{\mathrm{VM}}\right]$.

$${u}_{\mathrm{square}\text{\_}\mathrm{well}}={v}_{\mathrm{square}\text{\_}\mathrm{well}}\text{}\forall \text{}h{\sigma}_{\mathrm{square}\text{\_}\mathrm{well}}$$

$${u}_{\mathrm{ramp}}^{}={v}_{\mathrm{ramp}}\mathrm{max}\left(\frac{{\sigma}_{\mathrm{ramp}}-h}{{\sigma}_{\mathrm{ramp}}},0\right)$$

In absence of solid surfaces (red lines in Figure 3, θ = 180°), the profiles exhibit a characteristic sigmoid shape, whereas the corresponding surface tension becomes γ^{VM} ~73.0, ~12.0, and ~29.5 mJ/m^{2}, for HFD, SL, and SL-SGT, respectively [10,11]. In HFD, the isothermal compressibility has been set equal to κ_{T} = 1.43 GPa^{–1} [10,32], while the compressibility is roughly the same in the SL models (see Table S1 in Ref. [28]). The experimental surface tension of PE is 26.5–27.7 mJ/m^{2} [41], whereas the corresponding atomistic profile has a span of ~1 nm, [11]; thus, the SL-SGT model is more suitable for describing polymer–vacuum interfaces.

With increasing polymer–solid interactions (u_{s}) the profiles move closer to the solid surfaces and become more pronounced, especially when the HFD EoS is employed. Note that the profile with u_{square_well} for $\theta $= 45.3° is identical to the corresponding profiles in Refs. [10,32].

Even though the compressibilities of these models are very similar, the density profiles from SL and SL-SGT EoS are more expanded because they take into consideration the probable formation of gas phases. In addition, they are less pronounced due to the existence of a logarithmic term that suppresses large fluctuations of the density above unity. In SL-SGT, the profiles are almost identical for the same contact angles regardless of the functional form of the potential (square well versus ramp).

Besides analytic functional forms, RuSseL supports tabulated potentials as well. Even though such potentials might be more cumbersome to work with, they are more flexible, in the sense that, in many cases, they allow for the reproduction of profiles of arbitrary shape.

For example, one can optimize the tabulated potentials to reproduce target density profiles, φ_{target}, via iterative Boltzmann inversion:
with a being a relaxation parameter. This process can be envisioned as reverse engineering the self-consistent field; instead of trying to predict the density profiles for a given field, the optimizer attempts to find the field that reproduces the target density profiles.

$${u}_{\mathrm{s},i}\left(r\right)={u}_{\mathrm{s},i-1}\left(r\right)+a{k}_{\mathrm{B}}T\mathrm{ln}\left(\frac{{\phi}_{i-1}\left(r\right)}{{\phi}_{\mathrm{target}}\left(r\right)}\right)$$

Figure 4a depicts a density profile at a polyethylene–graphite interface at 450 K obtained from molecular dynamics simulations [33] and the optimized density profiles from RuSseL using the HFD and the ideal gas (IG) free energy density; in the latter, the intermolecular interactions among chains are turned off, while the chain segments interact explicitly with the solid walls. The corresponding tabulated potentials are shown in Figure 4c. According to Figure 4a, it is possible to reproduce the MD profiles exactly, given that the underlying EoS does not impose any particular constraints. For example, it is impossible to reproduce these profiles with SL EoS since the logarithmic term does not allow the density to exceed the characteristic SL density. Similarly, Figure 4b depicts an exotic target density profile where φ = 0 for h $\in $ [1, 2), φ = 1.1 for h $\in $ [2, 3), and φ = 1 everywhere else. As with the previous case, the target profile has been reproduced with HFD, SL, and IG EoS in the presence of the corresponding tabulated potentials shown in Figure 4d.

The first step for the mathematical representation of the problem is to define the Hamiltonian of the system as a functional of the chemical potential field w, the segment density field, ρ, and the partition functions of matrix, Q_{m}, and grafted chains, Q_{g} relative to the field-free case. The methodology can also address systems comprising exclusively chains grafted on a solid surface, that is, not in contact with polymer melt. This is achieved by switching to a canonical ensemble formalism [28]. In this work, we limited the presentation of the mathematical formulation to cases where polymer melt is present. When matrix chains are present in the system, the grand canonical ensemble is most suitable for the thermodynamic description and the Hamiltonian, H, is given by Equation (4).
where μ_{m} is the chemical potential per segment of matrix chains, u_{s} describes the polymer segment–solid interactions, and f is an excess free energy density functional, which is used to describe the nonbonded cohesive interactions among the polymer segments. V is the volume occupied by polymer in the system; n_{g} is the number of grafted chains, each grafted at position ${r}_{\mathrm{g},}{}_{{i}_{\mathrm{g}}}$; N_{m} and N_{g} are the lengths of matrix and grafted chains, in segments; Z_{m,free} and Z_{g,free} are partition functions of field-free matrix and grafted chains, respectively; β = 1/k_{B}T, where k_{B} is the Boltzmann constant and T the temperature of the polymer melt; and $A$_{m} and $A$_{g} are normalizing factors.

$$\begin{array}{l}H\left[\rho (r),\nabla \rho (r),w(r)\right]={\displaystyle \int \mathrm{d}r\left\{-iw(r)\rho (r)+f[\rho (r),\nabla \rho (r)]+{u}_{\mathrm{s}}(r)\rho (r)\right\}}\\ -\frac{1}{\beta}\mathrm{exp}(\beta {\mu}_{\mathrm{m}}{N}_{\mathrm{m}}){A}_{\mathrm{m}}{Z}_{\mathrm{m},\mathrm{free}}V{Q}_{\mathrm{m}}[iw(r)]-\frac{1}{\beta}{\displaystyle \sum _{{i}_{\mathrm{g}}=1}^{{n}_{\mathrm{g}}}\mathrm{ln}\left({A}_{\mathrm{g}}{Z}_{\mathrm{g},\mathrm{free}}{Q}_{\mathrm{g}}[{r}_{\mathrm{g},}{}_{{i}_{\mathrm{g}}};iw(r)]\right)}\end{array}$$

In order to reduce the complexities arising from the functional integration of the spatially varying fields w and ρ, it is plausible to neglect the fluctuations of the fields by invoking a saddle-point approximation (also known as “mean-field approximation”), and to maintain only the configuration of the fields with the maximum statistical weight in the expression for the grand potential in terms of the Hamiltonian of Equation (4) [42]. In doing so, the interactions between the polymer segments are replaced with the interaction of each individual segment with the mean-field, which is generated by the rest of the segments, and the segment density field is connected to the real field ${w}^{\prime}=iw$ (since w is purely imaginary) via Equation (5).

$${{w}^{\prime}}_{\mathrm{ifc}}(r)={w}^{\prime}(r)-{{w}^{\prime}}_{\mathrm{bulk}}={\frac{\partial f[\rho ,\nabla \rho ]}{\partial \rho}|}_{\rho =\rho (r)}-{\frac{\partial f[\rho ,\nabla \rho ]}{\partial \rho}|}_{\rho ={\rho}_{\mathrm{seg},\mathrm{bulk}}}-\nabla \cdot \frac{\partial f[\rho ,\nabla \rho ]}{\partial \nabla \rho}+{u}_{\mathrm{s}}(r)$$

Since the SCFT addresses systems of long polymer chains and high-density melts, it is quite common to represent the polymer chains with the Gaussian thread model. Daoulas et al. [32] have studied the accuracy of the Gaussian thread model in comparison to the worm-like chain model for the thermodynamic description of semi-flexible polymer melts at interfaces using SCFT. When chains are mathematically described as Gaussian threads, their conformations are governed by the Edwards diffusion Equation (6) [43,44].
where the solution, q_{c}, is the restricted partition function for chains of kind c [10,11,27,32,42]; R_{g,c} is the radius of gyration of a polymer chain of kind c; N is the variable spanning the contour of the chain; N_{c}, as in Equation (4), is the length of a chain of kind c in monomer segments; and ${{w}^{\prime}}_{\mathrm{ifc}}\left(r\right)$, as in Equation (5), is the self-consistent field at a certain point **r** of the domain minus the corresponding value of the field in the bulk polymer region [27]. The code offers the possibility to solve Equation (6) for any thermoplastic polymer, providing the appropriate value for the parameter of the “diffusion coefficient” ${R}_{\mathrm{g},c}{}^{2}/{N}_{c}$ and the parameters for the EoS used to calculate the excess free energy density, for both matrix and grafted chains; the difference between the two cases being that they require different initial conditions. In the case of matrix chains, the initial conditions across the domain are given by Equation (7):
Whereas, for the grafted chains, Equation (8) holds:
where ρ_{seg,bulk} is the density of polymer segments in the bulk region, σ_{g} is the areal density of grafted chains and ${h}_{\mathrm{g}}$ is the position of the grafting point of the grafted chain. The ratio of solid to grafting surface areas appearing in front of the right-hand side of eq 8 departs from unity only in nonplanar geometries. The value of the delta function is evaluated geometrically as the inverse of the interval length assigned to the grafting node of the 1D spatial mesh; for more details, see Section 2.1.2 of Ref. [27].

$$\frac{\partial}{\partial N}{q}_{c}(r,N)=\frac{{R}_{\mathrm{g},c}{}^{2}}{{N}_{c}}{\nabla}_{r}^{2}{q}_{c}(r,N)-\beta {{w}^{\prime}}_{\mathrm{ifc}}(r){q}_{c}(r,N)\text{}(c{=\mathrm{m},\text{}\mathrm{g}}^{-}{,\text{}\mathrm{g}}^{+})$$

$${q}_{\mathrm{m}}(r,0)=1$$

$${q}_{\mathrm{g}}({h}_{\mathrm{g}},0)=\frac{{S}_{\mathrm{solid}}}{{S}_{{h}_{\mathrm{g}}}}\frac{{\sigma}_{\mathrm{g}}{N}_{\mathrm{g}}}{{\rho}_{\mathrm{seg},\mathrm{bulk}}}\frac{\delta \left(h-{h}_{\mathrm{g}}\right)}{{q}_{\mathrm{m}}\left({h}_{\mathrm{g}},{N}_{\mathrm{g}}\right)}$$

In several works [10,18,32], the Crank–Nicholson contour discretization is implemented, which is a semi-implicit scheme (also known as “central differences”), in that the unknown solution, ${q}_{h}^{N}$, is expressed in our implementation by means of a central differences scheme, averaged between two successive contour points, N and N+ΔΝ, as shown in Equation (9).
while the first derivative of the solution, q, with respect to the contour variable N is given by Equation (10).
Therefore, the matrix form of the partial differential Equation (6) to be solved becomes as follows:
with $D=\frac{{R}_{\mathrm{G},c}{}^{2}\Delta N}{2{N}_{c}\Delta {h}^{2}}$ being the diffusion coefficient.

$$\frac{{\partial}^{2}q}{\partial {h}^{2}}=\frac{1}{2}\frac{{q}_{h+1}{}^{N+1}-2{q}_{h}{}^{N+1}+{q}_{h-1}{}^{N+1}}{\Delta {h}^{2}}+\frac{1}{2}\frac{{q}_{h+1}{}^{N}-2{q}_{h}{}^{N}+{q}_{h-1}{}^{N}}{\Delta {h}^{2}}$$

$$\frac{\partial {q}_{h}}{\partial N}=\frac{{q}_{h}{}^{N+1}-{q}_{h}{}^{N}}{\Delta N}$$

$$-D{q}_{h-1}{}^{N+1}+\left(1+2D+\frac{\Delta N\beta {w}_{\mathrm{ifc},h}}{2}\right){q}_{h}{}^{N+1}-D{q}_{h+1}{}^{N+1}=D{q}_{h-1}{}^{N}+\left(1-2D-\frac{\Delta N\beta {w}_{\mathrm{ifc},h}}{2}\right){q}_{h}{}^{N}+D{q}_{h+1}{}^{N}$$

In matrix-vector notation, Equation (11) can be written as presented in Equation (12).
where **I** is the identity matrix and we have defined the matrices **T** and **W**, which are shown below:

$$\left(I-DT+\frac{\Delta N}{2}W\right){q}^{N+1}=\left(I+DT-\frac{\Delta N}{2}W\right){q}^{N}$$

$$T=\left[\begin{array}{cccccc}-2& 1& & & & 1\\ 1& -2& 1& & & \\ & 1& -2& 1& & \\ & & \dots & \dots & \dots & \\ & & & 1& -2& 1\\ 1& & & & 1& -2\end{array}\right]$$

$$W=\left[\begin{array}{cccccc}w(0)& & & & & \\ & w(\Delta h)& & & & \\ & & w(2\Delta h)& & & \\ & & & \dots & & \\ & & & & w(L-\Delta h)& \\ & & & & & w(L)\end{array}\right]$$

Alternatively, Equation (12) can be written in terms of the stiffness matrix, $K=I-DT+\frac{\Delta N}{2}W$, and the vector denoting the right-hand side, $R=\left(I+DT-\frac{\Delta N}{2}W\right){q}^{N}$ (i.e., the solution vector at the previous contour-step weighted with $I+DT-\frac{\Delta N}{2}W$), as follows:
where ${a}_{i}=\u2013D$, ${b}_{i}=1+2D+\frac{\Delta N\beta {w}_{\mathrm{ifc},h}}{2}$ and ${c}_{i}=\u2013D$. Note that Equation (15) has been written in the general form, without imposing any boundary conditions; all nodes are equivalent and as a result, the domain is periodic.

$$K{q}^{N+1}=R\Rightarrow \left[\begin{array}{cccccc}{b}_{1}& {c}_{1}& & & & {a}_{1}\\ {a}_{2}& {b}_{2}& {c}_{2}& & & \\ & {a}_{3}& {b}_{3}& {c}_{3}& & \\ & & \dots & \dots & \dots & \\ & & & {a}_{n-1}& {b}_{n-1}& {c}_{n-1}\\ {c}_{n}& & & & {a}_{n}& {b}_{n}\end{array}\right]\left[\begin{array}{c}{q}_{1}^{N+1}\\ {q}_{2}^{N+1}\\ {q}_{3}^{N+1}\\ \dots \\ {q}_{n-1}^{N+1}\\ {q}_{n}^{N+1}\end{array}\right]=\left[\begin{array}{c}{R}_{1}^{}\\ {R}_{2}\\ {R}_{3}\\ \dots \\ {R}_{n-1}\\ {R}_{n}^{}\end{array}\right]$$

A more stable, but computationally more demanding, way to solve the time-dependent partial differential equation is to express the unknown solution, ${q}_{h}^{N}$, in terms of the next contour step, N+ΔN, according to Equation (16), presented below.

$$\frac{{\partial}^{2}q}{\partial {h}^{2}}=\frac{{q}_{h+1}{}^{N+1}-2{q}_{h}{}^{N+1}+{q}_{h-1}{}^{N+1}}{\Delta {h}^{2}}$$

In both cases, a linear system of equations needs to be solved to determine the solution at each contour step, but this implicit contour stepping method (also known as “backward differences”) allows for larger contour steps without reaching the numerical stability limits of the semi-implicit case. The first derivative of the chain contour is still given by Equation (10).

Adopting this discretization scheme, we end up with the following matrix-form of the partial differential equation to be solved:
where, again, the diffusion coefficient is given by the expression, $D=\frac{{R}_{\mathrm{G},c}{}^{2}\Delta N}{2{N}_{c}\Delta {h}^{2}}$. In matrix-vector notation, Equation (17) is written as:
where **I** is the identity matrix, and the matrices **T** and **W** are, again, those presented inEquation (14). In contrast to the semi-implicit scheme developed in the previous section, in the implicit one, the right-hand side is just the solution vector of the previous contour-step. As a result, eq 15 transforms as follows: ${a}_{i}=-2D$, ${b}_{i}=1+4D+\Delta N\beta {w}_{\mathrm{ifc},h}$, ${c}_{i}=-2D$, and $R={q}^{N}$.

$$-2D{q}_{h-1}{}^{N+1}+\left(1+4D+\Delta N\beta {w}_{\mathrm{ifc},h}\right){q}_{h}{}^{N+1}-2D{q}_{h+1}{}^{N+1}={q}_{h}{}^{N}$$

$$\left(I-2DT+\Delta NW\right){q}^{N+1}={q}^{N}$$

In formulating the matrices presented above, we have not taken into consideration the boundary conditions that occur from the physics of the problem. Given the single dimensionality of the systems addressed by the version of RuSseL presented in this work, the system is mathematically represented by a line, which is bounded by one point on the left and another one on the right. In most cases where aperiodic systems are addressed, at least one of the bounding points needs to be assigned a Dirichlet boundary condition (also known as an “essential” or “absorbing” boundary condition). The physical interpretation of this type of condition is that the polymer melt is in contact with a solid or gas surface, and thus, the polymer segments are not allowed to reach that surface. On the other point of the domain, we can assign either a Dirichlet boundary condition as well, or a Neumann boundary condition and, therefore, specify a certain value for the derivative of the solution rather than the solution itself. If the right-hand side point denotes the position where the bulk region starts and the system is symmetric, then the solution derivative is set equal to zero. In order to demonstrate these considerations, we present the linear system of equations to be solved, in the case where Dirichlet or Neumann boundary conditions are imposed on both boundary points. For simplicity, we impose the boundary conditions at the nodes lying at the edges of the domain; it is possible, however, to apply them anywhere across the domain, as we demonstrate in Section 5.3. and 5.4.

In a system with Dirichlet BC at i = 1 and i = n, eq 15 becomes the following:

$$\left[\begin{array}{cccccc}1& 0& & & & 0\\ {a}_{2}& {b}_{2}& {c}_{2}& & & \\ & {a}_{3}& {b}_{3}& {c}_{3}& & \\ & & \dots & \dots & \dots & \\ & & & {a}_{n-1}& {b}_{n-1}& {c}_{n-1}\\ 0& & & & 0& 1\end{array}\right]\left[\begin{array}{c}{q}_{1}\\ {q}_{2}\\ {q}_{3}\\ \dots \\ {q}_{n-1}\\ {q}_{n}\end{array}\right]=\left[\begin{array}{c}{R}_{1}^{\mathrm{DIR}}\\ {R}_{2}\\ {R}_{3}\\ \dots \\ {R}_{n-1}\\ {R}_{n}^{\mathrm{DIR}}\end{array}\right]$$

In practice, applying the Dirichlet BC to the i^{th} node entails the following substitutions: a_{i} = 0, b_{i} = 1, c_{i} = 0, and ${R}_{i}^{}={R}_{i}^{\mathrm{DIR}}$ (corresponding to a fixed q_{i} value).

If zero-derivative Neumann boundary conditions are imposed on the left (i = 1) and right (i = n) edges of the domain, then the matrix **T** needs to be modified, which, in turn, influences the final stiffness matrix and the right-hand side vector.

For the semi-implicit scheme, the matrix **T** is given by Equation (20):
Whereas, for the implicit scheme, it is given by Equation (21).

$$T=\left[\begin{array}{cccccc}-1& 1& & & & 0\\ 1& -2& 1& & & \\ & 1& -2& 1& & \\ & & \dots & \dots & \dots & \\ & & & 1& -2& 1\\ 0& & & & 1& -1\end{array}\right]$$

$$T=\left[\begin{array}{cccccc}-2& 2& & & & 0\\ 1& -2& 1& & & \\ & 1& -2& 1& & \\ & & \dots & \dots & \dots & \\ & & & 1& -2& 1\\ 0& & & & 2& -2\end{array}\right]$$

In the case of non-periodic systems (a_{1} = c_{n} = 0), where the stiffness matrix assumes a tridiagonal form, the linear system of equations is solved with the Thomas algorithm [45]. As in the conventional lower-upper (LU) decomposition algorithm, the solution of the tridiagonal system comprises three essential steps: decomposition, forward substitution, and backward substitution, which are presented below:

Decomposition:
Forward substitution:
Backward substitution:

$$\begin{array}{l}{a}_{i}\to \frac{{a}_{i}}{{b}_{i-1}}\\ {b}_{i}\to {b}_{i}-{a}_{i}{c}_{i-1}\text{},\text{}i=2,\dots ,n-1\end{array}$$

$${R}_{i}\to {R}_{i}-{a}_{i}{R}_{i-1},\text{}i=2,\dots ,n$$

$$\begin{array}{l}{q}_{n}^{N+1}\to \frac{{R}_{n}}{{b}_{n}}\\ {q}_{i}^{N+1}\to \frac{{R}_{i}-{c}_{i}{q}_{i+1}^{N+1}}{{b}_{i}},\text{}i=n-1,\dots ,1\end{array}$$

On the other hand, when considering periodic systems, the linear system of equations is solved with a more general (but computationally more expensive) solver, which is based on the traditional Gauss elimination method with pivoting functionality [46].

The code was written in Fortran95 and it can be compiled and run both on Windows and Linux systems. The code is modular, extensible, and the subroutines have been abstracted in a way that makes them reusable. The models and calculation parameters are entirely controlled through an input file, which is parsed during runtime and contains information concerning the polymer physics, the discretization of space and chain contour, and the convergence parameters.

During the initialization phase, the code attempts to read an input script by a dedicated subroutine, which parses the input file line-by-line and searches for special IDENTIFIERS that are related to specific user commands and operations. Whenever an IDENTIFIER is spotted, the program attempts to record the arguments(s) that precede it. In cases where an input line lacks any special IDENTIFIER or starts from “#” or “!”, it is skipped; thus, with very few exceptions, the order of the commands does not matter.

In addition, the parser has been equipped with an error-checking section that checks (i) whether the user has specified all the necessary inputs (e.g., domain size (L), temperature (T), geometry), (ii) the viability of individual user inputs (e.g., “are T or ρ_{mass,bulk} > 0?”, or “does a specific option exist?”), and (iii) whether the combinations of the inputs are viable (e.g., “is the grafting point located above the reflective wall?”). In cases where it finds a problem with the user inputs, the program issues a relevant error message and terminates.

The user has to initially specify the parameters of the domain, such as its size, geometry (planar, spherical), and boundary conditions (e.g., Dirichlet or Neumann), as well as other model-specific parameters, such as whether matrix or grafted chains are present in the system and the corresponding chain lengths, and thermodynamic conditions, such as the temperature and pressure. For spherical geometries the radius of the sphere, R_{NP}, is also required.

The mass of the polymer chains and the coarse-graining degree are set by the molecular weight of the monomer constituting the polymer chains, m_{monomer}, whereas the radius of gyration of the chains is set indirectly based on the following relation [44]:
with C_{∞} being the characteristic ratio and l_{C-C} being the chemical bond length between consecutive polymer segments provided by the user.

$$\frac{\langle {R}_{\mathrm{g},c}{}^{2}\rangle}{{N}_{c}}=\frac{1}{6}{C}_{\infty}{l}_{\mathrm{C}-\mathrm{C}}{}^{2}$$

If calculations regarding solid–polymer interfaces are to be performed, then the user must select the sides of the walls (low, high, or both) and the functional form of the potentials that describe these interactions. So far, the code supports the following potential styles: Vacuum (no potential), Hamaker potential (sphere–sphere, sphere–plate, plate–plate) [27,35,36,47], Square well, Ramp, Tabulated, Hybrid, and Custom. The Tabulated potential allows the user to input a potential table that includes the polymer–solid interactions as a function of the distance of a segment from the solid surface. The hybrid style allows the user to enable several functional forms at the same time. Finally, the custom style allows the user to implement a custom functional form by specifying the list of coefficients in the input file. The functional form, however, must be specified in compile time. For more details, see Supplementary Information.

The cohesion of the polymer is specified by assigning an appropriate EoS coupled with its necessary parameters. So far, the code offers the possibility to run SCFT calculations using the Helfand (HFD) [48] and the Sanchez–Lacombe (SL) [49] equations of state. In addition, the free energy densities can be combined with a square gradient term to more accurately address polymer–gas interfaces [11,50]. Nevertheless, the code has been written in a generic way, so that any other appropriate model can be inserted and used [27]. In general, our code is able to run SCFT calculations for any thermoplastic polymer melt, given that an appropriate EoS coupled with its necessary parameters is provided.

The stability of the iterative convergence scheme can be enhanced by adjusting several parameters such as the following:

- the contour length and spatial discretization (Δh, ΔΝ
_{c}) - the kind of discretization (uniform, nonuniform)
- the integration method (Simpson’s method, rectangle rule)
- the field relaxation parameter a
_{mix}, which depends on the isothermal compressibility, κ_{T}, of the bulk polymer and the length of the polymer chains [27] - the chain contour stepping method for the solution of the time-dependent PDE must be specified, that is, semi-implicit (also known as “Crank–Nicholson”) or implicit (see Section 3.2.1 and Section 3.2.2) [51].

The convergence criteria are based on the maximum number of iterations and the tolerance in the maximum error of the field, $\Delta {{w}^{\prime}}_{\mathrm{ifc}}^{\mathrm{tol}}$.

The flow chart in Figure 5 illustrates the structure and overall flow of the RuSseL code. The code can be split into two main sections: (i) the initialization sections, and (ii) the iterative convergence procedure.

During the initialization phase, the code parses the input file and checks for user errors. When the HFD EoS is selected, ρ_{mass,bulk} and κ_{T} are assigned the input values at the given temperature and pressure, whereas, for the SL-EoS, these are estimated automatically based on the Sanchez–Lacombe vapor–liquid equilibrium (SLVLE) of the fluid and the input characteristic density, temperature, and pressure (see SI in Ref. [11]). After allocating and initializing the required arrays, the program implements the spatial and contour length discretization based on the user inputs. Subsequently, the code sets the wall–segment distance based on the coordinates of the hard-sphere wall (either via the input values or based on the automatic recalibration) and then computes the surface area of each layer S(h) and the system volume, V, based on the specified geometry. The program then sets the kind and the coefficients of the polymer–wall interactions of each solid surface, and then initializes the self-consistent field, either to zero or by reading it from an input binary file (useful for restarting a calculation).

Subsequently, the program enters the iterative convergence procedure until the convergence criteria are satisfied. At the start of each cycle and for each kind of chain (c = m, g^{–}, g^{+}), the program computes the corresponding restricted partition function and reduced density as follows: (i) the partition function is initialized based on the boundary conditions and on Equation (7) or Equation (8) for matrix or grafted chains, respectively; (ii) the Edwards diffusion equation is evaluated for N up to N_{c}, subject to the field, ${{w}^{\prime}}_{\mathrm{ifc}}$, in planar or spherical geometries; (iii) the reduced density of the kind-c chains is calculated as follows:
with C being a convolution integral of the form:

$${\phi}_{c}(r)=C\left({q}_{c},{q}_{\mathrm{m}},{N}_{c},r\right)$$

$$C\left({q}_{\mathrm{a}},{q}_{\mathrm{b}},{N}_{c},r\right)=\frac{1}{{N}_{\mathrm{c}}}{\displaystyle \underset{0}{\overset{{N}_{\mathrm{c}}}{\int}}\mathrm{d}N}\text{}{q}_{\mathrm{a}}(r,N)\text{}{q}_{\mathrm{b}}(r,{N}_{\mathrm{c}}-N)$$

Next, the new field, ${{w}^{\prime}}_{\mathrm{ifc}}^{\mathrm{new}}$, is calculated from Equation (5) from the total reduced density, $\phi ={\phi}_{\mathrm{m}}+{\phi}_{{\mathrm{g}}^{-}}+{\phi}_{{\mathrm{g}}^{+}}$, and is then relaxed based on the following mixing rule:

$${{w}^{\prime}}_{\mathrm{ifc}}(r)\to \left(1-{a}_{\mathrm{mix}}\right){{w}^{\prime}}_{\mathrm{ifc}}(r)+{a}_{\mathrm{mix}}{{w}^{\prime}}_{\mathrm{ifc}}^{\mathrm{new}}(r)$$

At a frequency that is specified by the user, the program exports the thermodynamic properties and computes. Finally, in case the convergence conditions are satisfied ($\Delta {{w}^{\prime}}_{\mathrm{ifc}}^{\mathrm{tol}}<\mathrm{max}\left(\{\Vert {{w}^{\prime}}_{\mathrm{ifc}}^{\mathrm{new}}(r)-{{w}^{\prime}}_{\mathrm{ifc}}(r)\Vert ,\forall r\in R\}\right)\text{}\mathrm{or}\text{}{N}_{\mathrm{step}}{N}_{\mathrm{step}}^{\mathrm{max}}$), the program exports the outputs and finalizes, or else the cycle of the iterative procedure restarts.

During the calculation, RuSseL offers the possibility to compute several properties regarding the thermodynamics of the system and the configurations of matrix and grafted chains. The export frequency of each compute can be adjusted by the user in order to save computational time and hard disk space. Unless otherwise stated, the present section will illustrate the exported quantities from a perfectly wetted PS-Silica system (same as Figure 1) in a GMV geometry with ${\sigma}_{{\mathrm{g}}^{-}}$ = 0.4 nm^{–2}, ${N}_{{\mathrm{g}}^{-}}$ = 50, N_{m} = 100, and L = 10 nm.

To begin with, useful quantities are the total (φ_{total}) and partial reduced density profiles of the grafted (φ_{g}) and matrix (φ_{m}), such as those shown in Figure 1, Figure 2 and Figure 6 for the GMV geometry examined here. With φ known, the segmental or the mass density profiles can be retrieved by multiplying φ with ρ_{seg,bulk} or ρ_{mass,bulk} of the corresponding polymer, respectively.

The dimensions of the grafted chains can be quantified in terms of the root mean square brush thickness:
and in terms of ⟨h_{g,99%}⟩, which is defined here as the distance from the solid surface, which encloses 99% of the grafted chain segments:

$${\langle {h}_{\mathrm{g}}{}^{2}\rangle}^{1/2}={\left[\frac{{\displaystyle {\int}_{R}^{}\mathrm{d}h\text{}{h}^{2}{\rho}_{\mathrm{g}}(h)}}{{\displaystyle {\int}_{R}^{}\mathrm{d}h{\rho}_{\mathrm{g}}(h)}}\right]}^{1/2}$$

$$\int}_{{R}_{99\%}^{}}^{}\mathrm{d}h{\rho}_{\mathrm{g}}(h)}=0.99{N}_{g}{n}_{\mathrm{g$$

Our code allows for the decomposition of the density profiles into contributions of individual chain segments, such as chain ends and middle segments [11,33,52]. The contribution of the N^{th} segment to the corresponding density profile of the kind-c chain can be retrieved by the following equation:

$${\phi}_{c,N}(r)=\frac{1}{{N}_{c}}{q}_{c}(r,N)\text{}{q}_{\mathrm{m}}(r,{N}_{c}-N)$$

Setting N to 1 or N_{c} results in the density profile of the end segments, which can provide a very useful measure regarding the tendencies of the chain ends to segregate at the interfaces [11,27,33]. Setting N = N_{c}/2, on the other hand, results in the reduced density profiles of middle segments; comparisons among these two provide useful information regarding the overall shape of the chains [11].

Figure 7a depicts the segmental density profiles, φ_{c}_{,N}, for the chain ends (N = 1, N_{c}) and middle segments (N = N_{c}/2) of grafted and matrix chains. The first segment of the grafted chains (N = 1) corresponds to the grafting point and features a sharp peak at h = h_{g}. The middle segment and the last end segment of the grafted chains exhibit continuous density profiles, on the other hand, wherein the latter spreads further towards the matrix region. The density profiles of the matrix segments are suppressed in the vicinity of the grafted chains, since the latter reduce the available accessible space. Unlike grafted chains, the two end segments of the matrix chains are equivalent due to symmetry, hence their corresponding distributions are identical, ${\phi}_{\mathrm{m},1}={\phi}_{\mathrm{m},{N}_{\mathrm{m}}}$. The ends of the matrix chains feature more pronounced profiles near the interfaces as compared to the middle segments.

These tendencies of the chain segments to segregate at the interfaces can be better quantified in terms of the normalized segment distribution, which is calculated via Equation (32).

$${\tilde{\phi}}_{c,N}(r)=\frac{{N}_{c}{\phi}_{c,N}(r)}{{\phi}_{c}(r)}$$

In practice, ${\tilde{\phi}}_{c,N}(r)$ denotes the density of the Nth chain segments at a distance **r**, relative to the total segment density of the type-c chains; consequently, in a bulk phase, ${\tilde{\phi}}_{c,N}=1$.

According to Figure 7b, the segment density profiles of matrix chains become ${\tilde{\phi}}_{\mathrm{m},N}=1$, across the bulk region. Near the solid–polymer and the polymer–vacuum interfaces, the profiles of the end segments of the matrix chains are enhanced significantly by ${\tilde{\phi}}_{\mathrm{m},1}$~ 6 and ${\tilde{\phi}}_{\mathrm{m},{N}_{\mathrm{m}}}$~ 100, respectively, whereas the profiles of the middle segments are slightly less than bulk, ${\tilde{\phi}}_{\mathrm{m},{N}_{\mathrm{m}}/2}<1$. The corresponding profiles for the two ends of the grafted chains are highly asymmetric since ${\tilde{\phi}}_{\mathrm{g},1}$~210 and ${\tilde{\phi}}_{\mathrm{g},{N}_{\mathrm{g}}}~0.6$ at the grafting point, and ${\tilde{\phi}}_{\mathrm{g},1}$ = 0 and ${\tilde{\phi}}_{\mathrm{g},{N}_{\mathrm{g}}}~180$ at the edge of the free surface. The middle segments of the grafted chains exhibit suppressed profiles at the interfaces and over most of the matrix region. This effect is much more pronounced than for matrix chains; at the polymer–vacuum interface, ${\tilde{\phi}}_{\mathrm{g},{N}_{\mathrm{g}}/2}\ll 1$.

Besides the end and middle segments, RuSseL allows for the exporting of the profile of any other segment specified by the user. The contour plots in Figure 8 depict ${\phi}_{c,N}$ and ${\tilde{\phi}}_{c,N}$ for all the segments of grafted and matrix chains.

RuSseL has the option to decompose the density profiles into contributions of adsorbed and free segments based on segment surface distance criteria. Essentially, whenever a chain segment lies at a distance lower than ${h}_{\mathrm{ads}}^{\mp}$ from the left (–) or right (+) surface, it can be classified as adsorbed to that surface. This is a purely geometric distinction based on a critical distance from the solid surface, which is defined by the user through the input file. There are several approaches to setting the critical distance, depending on the specific application:

- Solid adsorption: h
_{ads}can be tuned based on the peaks of the density profile (e.g., in Ref. [37], h_{ads}was set equal to 0.6 nm, the distance between the first two peaks of the polyethylene/graphite density profile), or based on the strength of polymer–solid interactions (e.g., in Ref. [28], h_{ads}was set equal to 1.28 nm, where the PS-Silica interactions, as described by the Hamaker potential, become very weak). - Segregation at polymer–vacuum interfaces: In Ref. [11], h
_{ads}was set equal to a distance where the reduced density φ reaches 0.5. - Brush penetration: h
_{ads}can also be set to the span of the grafted brush, h_{g,99%}, in order to quantify the tendency of matrix chains to penetrate the brushes, or the tendencies of opposing brushes to penetrate each other.

Based on the distribution of adsorbed and free segments of a chain, it is possible to classify the chain into several states and sub-states, which are denoted in Table 1.

A chain that is comprised entirely of free segments is classified as free (f); otherwise, in case it includes adsorbed segments, it is treated as adsorbed (${\mathrm{a}}^{\mp}$). The adsorbed chains can be classified into fully (${\mathrm{a}}_{\mathrm{full}}^{\mp}$) and partially (${\mathrm{a}}_{\mathrm{part}}^{\mp}$) adsorbed chains in cases where they are entirely or partially adsorbed to a surface. Furthermore, the free segments of partially adsorbed chains can be classified into free loops (${\mathrm{a}}_{\mathrm{loop}\text{\_}\mathrm{f}}^{\mp}$) and tails (${\mathrm{a}}_{\mathrm{tail}\text{\_}\mathrm{f}}^{\mp}$) and the adsorbed ones as adsorbed loops (${\mathrm{a}}_{\mathrm{loop}\text{\_}\mathrm{a}}^{\mp}$) and tails (${\mathrm{a}}_{\mathrm{tail}\text{\_}\mathrm{a}}^{\mp}$); the latter are also classified into trains [8,11,32,37]. Finally, a chain that is partially adsorbed in both surfaces can be classified as a bridge, ${\mathrm{a}}_{\mathrm{bridge}}$. As can be imagined, the grafted chains cannot be classified to be free since, in most cases, the grafting points are located below h_{ads}; however, this procedure can unveil meaningful sub-states, such as fully and partially adsorbed grafted chains, as well as bridges.

Figure 9 includes some representative examples, whereas Figure 10 depicts the profiles of the aforementioned states from a perfectly wetted SMS system (see Table 1 in Ref. [27]) with ${h}_{\mathrm{ads}}^{-}={h}_{\mathrm{ads}}^{+}$ = 6 nm, and plate–plate distance, h_{ss}~ 16.8 nm (L = 16 nm). In this calculation, ${h}_{\mathrm{ads}}^{-}$ and ${h}_{\mathrm{ads}}^{+}$ have been chosen to be rather large, for illustrative purposes.

The first step of the classification procedure is to calculate the partition functions of the free (${q}_{c}^{\mathrm{f}}$), the fully adsorbed (${q}_{c}^{{\mathrm{a}}_{\mathrm{full}}^{-}},{q}_{c}^{{\mathrm{a}}_{\mathrm{full}}^{+}}$), and the non-adsorbed (${q}_{c}^{{!\mathrm{a}}^{-}},{q}_{c}^{{!\mathrm{a}}^{+}}$) chains with respect to the left (–) and the right (+) surface; for example, see Figure 10a. Note that a chain that is not adsorbed to a specific surface is not necessarily free since it might be adsorbed to the opposing surface.

To calculate each one of them, the Edwards diffusion equation is evaluated with the additional constraint that the Dirichlet boundary condition, q_{c}(h, N) = 0, is set to all the nodes with distance, h, within the ranges specified at the rightmost column of Table 1. In other words, the Edwards diffusion equation is evaluated in such a way that the chains are unable to access these ranges of the domain. Subsequently, the density profiles corresponding to these states and sub-states are calculated with the convolution integrals and the relations specified in Table 1. C symbolizes convolution of the restricted partition functions indicated; compare Equation (27).

Note that, ${\mathrm{a}}_{\mathrm{loop}}^{\mp}$ states corresponding to free or adsorbed segments are denoted as ${\mathrm{a}}_{\mathrm{loop}\text{\_}\mathrm{f}}^{\mp}$ and ${\mathrm{a}}_{\mathrm{loop}\text{\_}\mathrm{a}}^{\mp}$, respectively; hence, ${\mathrm{a}}_{\mathrm{loop}}^{\mp}={\mathrm{a}}_{\mathrm{loop}\text{\_}\mathrm{f}}^{\mp}+{\mathrm{a}}_{\mathrm{loop}\text{\_}\mathrm{a}}^{\mp}$. Moreover, note that the states have been defined in such a way that the following relations are satisfied:

$${\phi}_{c}^{{\mathrm{a}}^{\mp}}={\phi}_{c}^{{\mathrm{a}}_{\mathrm{full}}^{\mp}}+{\phi}_{c}^{{\mathrm{a}}_{\mathrm{part}}^{\mp}}$$

$${\phi}_{c}^{{\mathrm{a}}_{\mathrm{part}}^{\mp}}={\phi}_{c}^{{\mathrm{a}}_{\mathrm{loop}\text{\_}\mathrm{f}}^{\mp}}+{\phi}_{c}^{{\mathrm{a}}_{\mathrm{loop}\text{\_}\mathrm{a}}^{\mp}}+{\phi}_{c}^{{\mathrm{a}}_{\mathrm{tail}\text{\_}\mathrm{f}}^{\mp}}+{\phi}_{c}^{{\mathrm{a}}_{\mathrm{tail}\text{\_}\mathrm{a}}^{\mp}}$$

$${\phi}_{c}={\phi}_{c}^{{\mathrm{a}}^{-}}+{\phi}_{c}^{{\mathrm{a}}^{+}}+{\phi}_{c}^{\mathrm{f}}-{\phi}_{c}^{{\mathrm{a}}_{\mathrm{bridge}}}$$

Note that the bridges are a measure of the overlap between profiles of the adsorbed chains at the opposing surfaces; thus, in situations that the separation distance of these surfaces is large enough, φ_{bridge} → 0.

For both matrix and grafted chains, the number of chains per unit area can be determined, which is calculated through Equations (36) and (37). For the definition of this structural property, the reader is referred to previous works of authors [11,27,32,53].
${q}_{c,{h}_{0}}^{\mathrm{shape}}$ is the restricted partition function of all the chains that are unable to cross the surface at h_{0}, and it is calculated by evaluating the Edwards diffusion equation with the constraint that the Dirichlet boundary condition, q_{c}(h, N) = 0, is applied at the node at h_{0}. Subsequently, ${p}_{\mathrm{int},c}({h}_{0})$ is the probability that a chain will intersect this surface, and n_{ch,c}(h_{0}) corresponds to the number of type-c chains per unit area that pass through the surface at h_{0}. Consequently, near the grafting points, the chains/area equal the grafting density, as illustrated by the dashed line in Figure 11.

$${p}_{\mathrm{int},c}({h}_{0})=1-\frac{{\displaystyle \underset{R}{\overset{}{\int}}{q}_{c,{h}_{0}}^{\mathrm{shape}}(h,{N}_{c})\mathrm{d}h}}{{\displaystyle \underset{R}{\overset{}{\int}}{q}_{c}(h,{N}_{c})\mathrm{d}h}}$$

$${n}_{\mathrm{ch},c}({h}_{0})={p}_{\mathrm{int},c}({h}_{0})\frac{1}{{S}_{{h}_{0}}}\frac{1}{{N}_{c}}{\displaystyle \underset{R}{\overset{}{\int}}{\rho}_{c}(h)\mathrm{d}h}$$

The total free energy (grand potential) of the system is split into the following contributions:

$$\Delta \Omega =\Omega -{\Omega}_{\mathrm{bulk}}-{A}_{\mathrm{g},\mathrm{bulk}}=\Delta {\Omega}_{\mathrm{coh}}+\Delta {\Omega}_{\mathrm{field}}+\Delta {\Omega}_{\mathrm{m}}+\Delta {A}_{\mathrm{g}}+{U}_{\mathrm{s}}+{U}_{\mathrm{ss}}$$

Note that the total free energy has been defined with respect to a bulk melt of matrix chains with the same length (N_{m}) and chemical constitution, and at the same temperature and pressure, and with respect to ${n}_{{\mathrm{g}}^{\mp}}$ end-pinned unperturbed chains of length ${N}_{{\mathrm{g}}^{\mp}}$ exposed to bulk melt. In cases where the matrix chains are absent, Ω_{bulk} is zero. The cohesive interactions among the polymer segments (ΔΩ_{coh}) are described by Equation (39), with $f\left[\rho (r),\nabla \rho (r)\right]$ being a free energy density term from an equation of state specified by the user. The interaction between the (total) segment density field, ρ(**r**), and the self-consistent field, ${{w}^{\prime}}_{\mathrm{ifc}}\left(r\right)$, is given by Equation (40). The total solid–polymer interactions (U_{s})—when present—are given by Equation (41); whereas, their functional form, u_{s}(**r**), is specified by the user. The free energies associated with the conformational entropy of matrix and grafted chains (if they exist) are given by Equations (42) and (43), respectively.
where β = 1/k_{B}T. In systems with opposing solid surfaces, we introduce an additional enthalpic term that describes the solid–solid interactions. The functional form of this term depends on the potential of the solid.

$$\Delta {\Omega}_{\mathrm{coh}}={\displaystyle \underset{R}{\int}\mathrm{d}r}\left\{f\left[\rho (r),\nabla \rho (r)\right]-f\left[{\rho}_{\mathrm{seg},\mathrm{bulk}},0\right]\right\}$$

$$\Delta {\Omega}_{\mathrm{field}}=-{\displaystyle \underset{R}{\int}\mathrm{d}r}\left\{\rho (r){w}^{\prime}(r)-{\rho}_{\mathrm{seg},\mathrm{bulk}}{{w}^{\prime}}_{\mathrm{bulk}}(r)\right\}$$

$${U}_{\mathrm{s}}={\displaystyle \underset{R}{\int}\mathrm{d}r}\left\{\rho (r){u}_{\mathrm{s}}(r)\right\}$$

$$\Delta {\Omega}_{\mathrm{m}}=-\frac{{\rho}_{\mathrm{seg},\mathrm{bulk}}V}{\beta {N}_{\mathrm{m}}}\left({Q}_{\mathrm{m}}\left[{w}^{\prime}-{{w}^{\prime}}_{\mathrm{bulk}}\right]-1\right)$$

$$\Delta {A}_{\mathrm{g}}=-\frac{1}{\beta}{\displaystyle \sum _{{i}_{\mathrm{g}}=1}^{{n}_{\mathrm{g}}}\mathrm{ln}{Q}_{\mathrm{g}}\left[{r}_{\mathrm{g},{i}_{\mathrm{g}}};{w}^{\prime}-{{w}^{\prime}}_{\mathrm{bulk}}\right]}-\frac{1}{\beta}{\displaystyle \sum _{{i}_{\mathrm{g}}=1}^{{n}_{\mathrm{g}}}\mathrm{ln}\frac{{r}_{\mathrm{ref},q=0}^{}}{{r}_{\mathrm{g},{i}_{\mathrm{g}},q=0}^{}}}$$

An in-house code has been developed to perform calculations on one-dimensional domains based on Self-Consistent Field Theory and considering compressible polymer melts of Gaussian-thread chains. The code solves the Edwards diffusion equation using a Finite Differences scheme and it addresses heterogeneous polymer systems comprising homopolymer polymer melts in contact with gas or solid surfaces. In the latter case, the code offers the possibility to work with different potentials for the description of solid–polymer interactions, namely Hamaker, square-well, and user-specified potentials (either tabulated or analytic). There is also the option to choose from among the Helfand and the Sanchez–Lacombe equation of state, with or without square-gradient correction, for the description of the nonbonded interactions between the polymer segments. Building upon the content of our previous publications [11,27,28], we present the way that the code computes several thermodynamic and structural properties of the system.

The code has already been used by the authors in a number of publications regarding various systems, and this article, along with the accompanying manual, was written to serve as a technical reference for anyone who desires to reproduce the corresponding results. In this article, we have demonstrated a part of the capabilities of the code and validated our findings with results reported in the literature. There is plenty of room for future extensions of the code, for example, block copolymer systems, polydisperse polymer melts, additional models for nonbonded interactions (e.g., SAFT [54]), additional models for solid–polymer interactions, and investigation of more efficient solvers for the solution of the Edwards partial differential equations. The authors have also developed the three-dimensional analogue of RuSseL, implementing the Finite Element Method [25], which is also going to be published in the near future. The development of both versions of the code is an ongoing process and the authors are always trying to stay tuned and incorporate more computational tools and capabilities into their software, aiming to contribute to the area of polymer nanocomposites and shed light on several topics that still remain elusive to the community.

The source code of RuSseL and its manual are available online at https://0-www-mdpi-com.brum.beds.ac.uk/article/10.3390/computation9050057/s1.

Conceptualization, C.J.R., A.P.S. and D.N.T.; methodology, C.J.R., A.P.S., A.T.L. and D.N.T.; software, C.J.R., A.P.S. and A.T.L.; validation, C.J.R. and A.P.S.; formal analysis, C.J.R. and A.P.S.; investigation, C.J.R. and A.P.S.; resources, D.N.T.; data curation, C.J.R. and A.P.S.; writing—original draft preparation, C.J.R. and A.P.S.; writing—review and editing, D.N.T., C.J.R. and A.P.S.; visualization, A.P.S.; supervision, D.N.T.; project administration, D.N.T.; funding acquisition, D.N.T. All authors have read and agreed to the published version of the manuscript.

Financial support by the Hellenic Foundation for Research and Innovation (H.F.R.I), project number 1263, under the “First Call for H.F.R.I Research Projects to support Faculty members and Researchers and the procurement of high-cost equipment grant” is gratefully acknowledged. A.P.S. thanks SOLVAY, Paris, France, for financial support through the project “Multiscale Modeling for the Design of Antifouling Copolymers”. C.J.R. gratefully acknowledges financial support from the “Special Account for Research Funding” of the National Technical University of Athens and also through the project MuSiComPS, Grant Agreement 10062/19, sponsored by the LIMMAT Foundation, Switzerland. Part of this work was supported by computational time granted from the Greek Research and Technology Network (GRNET) in the National HPC facility—ARIS—under project ID pa201203 (MuSiPolI).

Not applicable.

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

- Yan, L.-T.; Xie, X.-M. Computational modeling and simulation of nanoparticle self-assembly in polymeric systems: Structures, properties and external field effects. Prog. Polym. Sci.
**2013**, 38, 369–405. [Google Scholar] [CrossRef] - Zeng, Q.H.; Yu, A.B.; Lu, G.Q. Multiscale modeling and simulation of polymer nanocomposites. Prog. Polym. Sci.
**2008**, 33, 191–269. [Google Scholar] [CrossRef] - Chen, S.; Doolen, G.D. Lattice Boltzmann Method for Fluid Flows. Annu. Rev. Fluid Mech.
**1998**, 30, 329–364. [Google Scholar] [CrossRef] - Español, P.; Warren, P.B. Perspective: Dissipative particle dynamics. J. Chem. Phys.
**2017**, 146, 150901. [Google Scholar] [CrossRef] - Bore, S.L.; Kolli, H.B.; De Nicola, A.; Byshkin, M.; Kawakatsu, T.; Milano, G.; Cascella, M. Hybrid particle-field molecular dynamics under constant pressure. J. Chem. Phys.
**2020**, 152. [Google Scholar] [CrossRef] [PubMed] - Caputo, S.; Hristov, V.; De Nicola, A.; Herbst, H.; Pizzirusso, A.; Donati, G.; Munaò, G.; Albunia, A.R.; Milano, G. Efficient Hybrid Particle-Field Coarse-Grained Model of Polymer Filler Interactions: Multiscale Hierarchical Structure of Carbon Black Particles in Contact with Polyethylene. J. Chem. Theory Comput.
**2021**, 17, 1755–1770. [Google Scholar] [CrossRef] - Sgouros, A.P.; Vogiatzis, G.G.; Megariotis, G.; Tzoumanekas, C.; Theodorou, D.N. Multiscale Simulations of Graphite-Capped Polyethylene Melts: Brownian Dynamics/Kinetic Monte Carlo Compared to Atomistic Calculations and Experiment. Macromolecules
**2019**, 52, 7503–7523. [Google Scholar] [CrossRef] - Sgouros, A.P.; Lakkas, A.T.; Megariotis, G.; Theodorou, D.N. Mesoscopic Simulations of Free Surfaces of Molten Polyethylene: Brownian Dynamics/kinetic Monte Carlo Coupled with Square Gradient Theory and Compared to Atomistic Calculations and Experiment. Macromolecules
**2018**, 51, 9798–9815. [Google Scholar] [CrossRef] - Lo Verso, F.; Egorov, S.A.; Milchev, A.; Binder, K. Spherical polymer brushes under good solvent conditions: Molecular dynamics results compared to density functional theory. J. Chem. Phys.
**2010**, 133, 184901. [Google Scholar] [CrossRef] [PubMed] - Theodorou, D.N.; Vogiatzis, G.G.; Kritikos, G. Self-consistent-field study of adsorption and desorption kinetics of polyethylene melts on graphite and comparison with atomistic simulations. Macromolecules
**2014**, 47, 6964–6981. [Google Scholar] [CrossRef] - Lakkas, A.T.; Sgouros, A.P.; Theodorou, D.N. Self-Consistent Field Theory Coupled with Square Gradient Theory of Free Surfaces of Molten Polymers and Compared to Atomistic Simulations and Experiment. Macromolecules
**2019**, 52, 5337–5356. [Google Scholar] [CrossRef] - Müller, M. Phase diagram of a mixed polymer brush. Phys. Rev. E Stat. Phys. Plasmas Fluids Relat. Interdiscip. Top.
**2002**, 65, 1–4. [Google Scholar] [CrossRef] - Trombly, D.M.; Ganesan, V. Curvature effects upon interactions of polymer-grafted nanoparticles in chemically identical polymer matrices. J. Chem. Phys.
**2010**, 133, 154904. [Google Scholar] [CrossRef] [PubMed] - Müller, M.; Schmid, F. Incorporating fluctuations and dynamics in self-consistent field theories for polymer blends. Adv. Polym. Sci.
**2005**, 185, 1–58. [Google Scholar] - Ouaknin, G.; Laachi, N.; Delaney, K.; Fredrickson, G.H.; Gibou, F. Self-consistent field theory simulations of polymers on arbitrary domains. J. Comput. Phys.
**2016**, 327, 168–185. [Google Scholar] [CrossRef] - Arora, A.; Qin, J.; Morse, D.C.; Delaney, K.T.; Fredrickson, G.H.; Bates, F.S.; Dorfman, K.D. Broadly Accessible Self-Consistent Field Theory for Block Polymer Materials Discovery. Macromolecules
**2016**, 49, 4675–4690. [Google Scholar] [CrossRef] - Rasmussen, K.; Kalosakas, G. Improved numerical algorithm for exploring block copolymer mesophases. J. Polym. Sci. Part B Polym. Phys.
**2002**, 40, 1777–1783. [Google Scholar] [CrossRef] - Drolet, F.; Fredrickson, G.H. Combinatorial screening of complex block copolymer assembly with self-consistent field theory. Phys. Rev. Lett.
**1999**, 83, 4317–4320. [Google Scholar] [CrossRef] - Drolet, F.; Fredrickson, G.H. Optimizing chain bridging in complex block copolymers. Macromolecules
**2001**, 34, 5317–5324. [Google Scholar] [CrossRef] - Kim, J.U.; Matsen, M.W. Finite-stretching corrections to the Milner-Witten-Cates theory for polymer brushes. Eur. Phys. J. E
**2007**, 23, 135–144. [Google Scholar] [CrossRef] - Vigil, D.L.; García-Cervera, C.J.; Delaney, K.T.; Fredrickson, G.H. Linear Scaling Self-Consistent Field Theory with Spectral Contour Accuracy. ACS Macro Lett.
**2019**, 8, 1402–1406. [Google Scholar] [CrossRef] - Ackerman, D.M.; Delaney, K.; Fredrickson, G.H.; Ganapathysubramanian, B. A finite element approach to self-consistent field theory calculations of multiblock polymers. J. Comput. Phys.
**2017**, 331, 280–296. [Google Scholar] [CrossRef] - Huttom, D.V. Fundamentals of Finite Element Analysis; McGraw Hill: New York, NY, USA, 2004; ISBN 0-07-112231-1. [Google Scholar]
- Smith, I.M.; Griffiths, D.V.; Margetts, L. Programming the Finite Element Method, 5th ed.; Wiley: West Sussex, UK, 2014; ISBN 9781119973348. [Google Scholar]
- Revelas, C.J.; Sgouros, A.P.; Lakkas, A.T.; Theodorou, D.N. A Three-Dimensional Finite Element Methodology for Addressing Heterogeneous Polymer Systems with Simulations Based on Self-Consistent Field Theory. In Proceedings of the International Conference of Computational Methods In Science and Engineering 2020 (ICCMSE 2020), Heraklion, Crete, Greece, 27 April–3 May 2020; p. 130002. [Google Scholar]
- Daoulas, K.C.; Müller, M. Exploring thermodynamic stability of the stalk fusion-intermediate with three-dimensional self-consistent field theory calculations. Soft Matter
**2013**, 9, 4097–4102. [Google Scholar] [CrossRef] - Lakkas, A.T.; Sgouros, A.P.; Revelas, C.J.; Theodorou, D.N. Structure and Thermodynamics of Grafted Silica/Polystyrene Nanocomposites Investigated Through Self-Consistent Field Theory. Soft Matter
**2021**, 17, 4077–4097. [Google Scholar] [CrossRef] - Sgouros, A.P.; Revelas, C.J.; Lakkas, A.T.; Theodorou, D.N. Potential of Mean Force between Bare or Grafted Silica/Polystyrene Surfaces from Self-Consistent Field Theory. Polymers (Basel)
**2021**, 13, 1197. [Google Scholar] [CrossRef] - Cheong, G.K.; Chawla, A.; Morse, D.C.; Dorfman, K.D. Open-source code for self-consistent field theory calculations of block polymer phase behavior on graphics processing units. Eur. Phys. J. E
**2020**, 43, 15. [Google Scholar] [CrossRef] - Price-Whelan, A.M. SCF 1991. Available online: https://github.com/adrn/scf_fortran (accessed on 10 May 2021).
- Bilchak, C.R.; Buenning, E.; Asai, M.; Zhang, K.; Durning, C.J.; Kumar, S.K.; Huang, Y.; Benicewicz, B.C.; Gidley, D.W.; Cheng, S.; et al. Polymer-Grafted Nanoparticle Membranes with Controllable Free Volume. Macromolecules
**2017**, 50, 7111–7120. [Google Scholar] [CrossRef] - Daoulas, K.C.; Theodorou, D.N.; Harmandaris, V.A.; Karayiannis, N.C.; Mavrantzas, V.G. Self-consistent-field study of compressible semiflexible melts adsorbed on a solid substrate and comparison with atomistic simulations. Macromolecules
**2005**, 38, 7134–7149. [Google Scholar] [CrossRef] - Sgouros, A.P.; Vogiatzis, G.G.; Kritikos, G.; Boziki, A.; Nikolakopoulou, A.; Liveris, D.; Theodorou, D.N. Molecular Simulations of Free and Graphite Capped Polyethylene Films: Estimation of the Interfacial Free Energies. Macromolecules
**2017**, 50, 8827–8844. [Google Scholar] [CrossRef] - Theodorou, D.N. Variable-Density Model of Polymer Melt/Solid Interfaces: Structure, Adhesion Tension, and Surface Forces. Macromolecules
**1989**, 22, 4589–4597. [Google Scholar] [CrossRef] - Vogiatzis, G.G.; Theodorou, D.N. Structure of polymer layers grafted to nanoparticles in silica-polystyrene nanocomposites. Macromolecules
**2013**, 46, 4670–4683. [Google Scholar] [CrossRef] - Hamaker, H.C. The London—van der Waals attraction between spherical particles. Physica
**1937**, 4, 1058–1072. [Google Scholar] [CrossRef] - Daoulas, K.C.; Harmandaris, V.A.; Mavrantzas, V.G. Detailed Atomistic Simulation of a Polymer/Solid Interface: Structure, Density and Conformation of a Thin Film of Polyethylene Melt Adsorbed on Graphite. Macromolecules
**2005**, 38, 5780–5795. [Google Scholar] [CrossRef] - Magda, J.J.; Tirrell, M.; Davis, H.T. Molecular dynamics of narrow, liquid-filled pores. J. Chem. Phys.
**1985**, 83, 1888–1901. [Google Scholar] [CrossRef] - Bilchak, C.R.; Jhalaria, M.; Huang, Y.; Abbas, Z.; Midya, J.; Benedetti, F.M.; Parisi, D.; Egger, W.; Dickmann, M.; Minelli, M.; et al. Tuning Selectivities in Gas Separation Membranes Based on Polymer-Grafted Nanoparticles. ACS Nano
**2020**, 14, 17174–17183. [Google Scholar] [CrossRef] - Mansfield, K.F.; Theodorou, D.N. Atomistic Simulation of a Glassy Polymer/Graphite Interface. Macromolecules
**1991**, 24, 4295–4309. [Google Scholar] [CrossRef] - Hong, K.M.; Noolandi, J. Conformational Entropy Effects in a Compressible Lattice Fluid Theory of Polymers. Macromolecules
**1981**, 14, 1229–1234. [Google Scholar] [CrossRef] - Fredrickson, G.H. The Equilibrium Theory of Inhomogeneous Polymers; International Series of Monographs on Physics; Oxford University Press: Oxford, UK, 2006. [Google Scholar]
- Doi, M.; Edwards, S.F. The Theory of Polymer Dynamics; Clarendon Press: Oxford, UK, 1986. [Google Scholar]
- Theodorou, D.N. Polymers at Surfaces and Interfaces. Comput. Simul. Surf. Interfaces
**2003**, 329–419. [Google Scholar] [CrossRef] - Higham, N.J. Accuracy and Stability of Numerical Algorithms; Society for Industrial and Applied Mathematics: University City, PA, USA, 2002; ISBN 978-0-89871-521-7. [Google Scholar]
- Datta, B.N. Numerical Linear Algebra and Applications, 2nd ed.; Society for Industrial and Applied Mathematics: Philadelphia, PA, USA, 2010; ISBN 0898716853. [Google Scholar]
- Vogiatzis, G.G.; Voyiatzis, E.; Theodorou, D.N. Monte Carlo simulations of a coarse grained model for an athermal all-polystyrene nanocomposite system. Eur. Polym. J.
**2011**, 47, 699–712. [Google Scholar] [CrossRef] - Helfand, E.; Tagami, Y. Theory of the interface between immiscible polymers. II. J. Chem. Phys.
**1972**, 56, 3592–3601. [Google Scholar] [CrossRef] - Poser, C.I.; Sanchez, I.C. Surface tension theory of pure liquids and polymer melts. J. Colloid Interface Sci.
**1979**, 69, 539–548. [Google Scholar] [CrossRef] - Lin, H.; Duan, Y.Y.; Min, Q. Gradient theory modeling of surface tension for pure fluids and binary mixtures. Fluid Phase Equilib.
**2007**, 254, 75–90. [Google Scholar] [CrossRef] - Press, W.H.; Teukolsky, S.A.; Vetterling, W.T.; Flannery, B.P. Numerical Recipes in C—The Art of Scientific Computing; Cambridge University Press: New York, NY, USA, 1992; ISBN 0521431085. [Google Scholar]
- Wu, D.T.; Fredrickson, G.H.; Carton, J.-P.; Ajdari, A.; Leibler, L. Distribution of chain ends at the surface of a polymer melt: Compensation effects and surface tension. J. Polym. Sci. Part B Polym. Phys.
**1995**, 33, 2373–2389. [Google Scholar] [CrossRef] - Theodorou, D.N. Lattice Models for Bulk Polymers at Interfaces. Macromolecules
**1988**, 21, 1391–1400. [Google Scholar] [CrossRef] - Chapman, W.G.; Gubbins, K.E.; Jackson, G.; Radosz, M. SAFT: Equation-of-state solution model for associating fluids. Fluid Phase Equilib.
**1989**, 52, 31–38. [Google Scholar] [CrossRef]

State | Symbol(α) | Reduced SegmentDensity | Restricted Part. Function | DirichletNodes, q_{c}(h, N) = 0 |
---|---|---|---|---|

free | f | ${\phi}_{c}^{\mathrm{f}}=C\left({q}_{c}^{\mathrm{f}},{q}_{c}^{\mathrm{f}},{N}_{c},h\right)$ | ${q}_{c}^{\mathrm{f}}$ | $h\in [0,{h}_{\mathrm{ads}}^{-}]\cup [{h}_{\mathrm{ads}}^{+},L]$ |

adsorbed fully | ${\mathrm{a}}_{\mathrm{full}}^{\mp}$ | ${\phi}_{c}^{{\mathrm{a}}_{\mathrm{full}}^{\mp}}=C\left({q}_{c}^{{\mathrm{a}}_{\mathrm{full}}^{\mp}},{q}_{c}^{{\mathrm{a}}_{\mathrm{full}}^{\mp}},{N}_{c},h\right)$ | ${q}_{c}^{{\mathrm{a}}_{\mathrm{full}}^{\mp}}$ | $h\in [{h}_{\mathrm{ads}}^{-},L]\text{}\mathrm{or}\text{}h\in [0,{h}_{\mathrm{ads}}^{+}]$ |

not adsorbed | ${!\mathrm{a}}^{\mp}$ | ${\phi}_{c}^{{!\mathrm{a}}^{\mp}}=C\left({q}_{c}^{{!\mathrm{a}}^{\mp}},{q}_{c}^{{!\mathrm{a}}^{\mp}},{N}_{c},h\right)$ | ${q}_{c}^{{!\mathrm{a}}^{\mp}}$ | $h\in [0,{h}_{\mathrm{ads}}^{-}]\text{}\mathrm{or}\text{}h\in [{h}_{\mathrm{ads}}^{+},L]$ |

adsorbed | ${\mathrm{a}}^{\mp}$ | ${\phi}_{c}^{{\mathrm{a}}^{\mp}}={\phi}_{c}-{\phi}_{c}^{{!\mathrm{a}}^{\mp}}$ | - | - |

adsorbed partially | ${\mathrm{a}}_{\mathrm{part}}^{\mp}$ | ${\phi}_{c}^{{\mathrm{a}}_{\mathrm{part}}^{\mp}}={\phi}_{c}^{{\mathrm{a}}_{}^{\mp}}-{\phi}_{c}^{{\mathrm{a}}_{\mathrm{full}}^{\mp}}$ | - | - |

loops | ${\mathrm{a}}_{\mathrm{loop}}^{\mp}$ | ${\phi}_{c}^{{\mathrm{a}}_{\mathrm{loop}}^{\mp}}=C\left({q}_{c}^{{\mathrm{a}}_{\mathrm{loop}}^{\mp}},{q}_{c}^{{\mathrm{a}}_{\mathrm{loop}}^{\mp}},{N}_{c},h\right)$ | ${q}_{c}^{{\mathrm{a}}_{\mathrm{loop}}^{\mp}}={q}_{c}-{q}_{c}^{{!\mathrm{a}}^{\mp}}-{q}_{c}^{{\mathrm{a}}_{\mathrm{full}}^{\mp}}$ | - |

tails outside the adsorbed region | ${\mathrm{a}}_{\mathrm{tail}\text{\_}\mathrm{f}}^{\mp}$ | ${\phi}_{c}^{{\mathrm{a}}_{\mathrm{tail}\text{\_}\mathrm{f}}^{\mp}}=2C\left({q}_{c}^{{\mathrm{a}}_{\mathrm{loop}}^{\mp}},{q}_{c}^{{!\mathrm{a}}_{}^{\mp}},{N}_{c},h\right)$ | - | - |

tails inside the adsorbed region | ${\mathrm{a}}_{\mathrm{tail}\text{\_}\mathrm{a}}^{\mp}$ | ${\phi}_{c}^{{\mathrm{a}}_{\mathrm{tail}\text{\_}\mathrm{a}}^{\mp}}=2C\left({q}_{c}^{{\mathrm{a}}_{\mathrm{loop}}^{\mp}},{q}_{c}^{{\mathrm{a}}_{\mathrm{full}}^{\mp}},{N}_{c},h\right)$ | - | - |

bridges | ${\mathrm{a}}_{\mathrm{bridge}}$ | ${\phi}_{c}^{{\mathrm{a}}_{\mathrm{bridge}}}={\phi}_{c}^{{\mathrm{a}}^{-}}+{\phi}_{c}^{{\mathrm{a}}^{+}}+{\phi}_{c}^{\mathrm{f}}-{\phi}_{c}$ | - | - |

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. |

© 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).