Description of interfaces of fluid-tethered chains: advances in density functional theories and off-lattice computer simulations

Many objects of nanoscopic dimensions involve fluid-tethered chain interfaces. These systems are of interest for basic science and for several applications, in particular for design of nanodevices for specific purposes. We review recent developments of theoretical methods in this area of research and in particular of density functional (DF) approaches, which provide important insights into microscopic properties of such interfaces. The theories permit to describe the dependence of adsorption, wettability, solvation forces and electric interfacial phenomena on thermodynamic states and on characteristics of tethered chains. Computer simulations for the problems in question are overviewed as well. Theoretical results are discussed in relation to simulation results and to some experimental observations.


Introduction
Systems involving interfaces formed by polymers tethered onto a solid surface and a fluid contacting them are widespread in nature and are the subject of many studies. One of the reasons is that the grafted layers of chains find several practical applications, for example in efficient colloid stabilization, protection and prevention of surfaces from undesirable adsorptions, manufacturing of functional surfaces and nanoscale fluidic channels, and many others. Interesting properties of such systems have also stimulated the development of several methods of synthesis for preparation of assemblies of end-grafted polymers. However, up to now, the development of new preparation methods is mainly based on empirical knowledge or intuition.
From the very beginning, experimental studies of tethered chains are accompanied by the development of theoretical methods of their description. They involve thermodynamic and structural properties of systems containing uncharged and charged species under the effect of various external factors including confinement at nanoscale as well.
The purpose of this article is to present a brief review of the statistical mechanical DF methods and offlattice computer simulations of classical fluids in contact with the walls modified by tethered polymers, to discuss some recent achievements in this area of research, as well as to provide some prospects for further investigations. After introductory remarks, in the next section (section 2), we briefly outline different types of theoretical approaches used for the description of tethered chains-fluid interfaces. We begin with the scaling theories, but our emphasis is on the theories that are based on the microscopic description and the use of the intermolecular potentials as the input, i.e., on the integral equations, self-consistent field and density functional theories. Particular attention is paid to the DF approaches. In section 3 we discuss one version of the density functional theory (DFT) that has been successful in describing several phenomena occurring within tethered chains-fluid interfaces. Some illustrative examples are given in section 4. Section 5 extends the approach of section 4 to the case of systems having electrostatic interactions. We discuss the theory and then, in section 6, present some applications concerning structural, thermodynamic and electric properties of the systems with tethered chains. Mesoscopic simulations by the dissipative particle dynamics method and the (semi-) atomistic simulations using the Monte Carlo (MC) and Molecular Dynamics (MD) techniques are reviewed in sections 7.1 and 7.2, respectively. The manuscript ends with a short summary in section 8.

Overview of theoretical methods
Essential features of the systems involving end-tethered polymers onto a solid and in contact with a fluid can be captured by modelling at atomistic or semi-atomistic (at groups of atoms) scale by using contemporary theories for interfaces between a solid and a complex fluid. A link between fundamentals and applications can be reached by combining theoretical (classical and quantum) tools, intensive computer simulations and laboratory experiments. Current studies at a simple level of modelling put focus on the development of a theory for various phenomena occurring in the systems of our interest. Computer simulations provide data to confront theoretical predictions and simulation results. The most important, they serve as predictive tools to guide elaboration of interfaces with desired functional properties.
Theoretical tools to investigate systems with tethered chains are similar to the methods designed to study nonuniform polymeric systems. They can be classified into several groups by using different criteria. One can divide them according to the way in which translational degrees of freedom are considered. Namely, either lattice or off-lattice models can be distinguished. On the other hand, the classification can be made according to theoretical approaches used to describe a model. One can distinguish the theories based on the application of scaling arguments [1,2], integral equations approaches, self-consistent field (SCF) as well as the DF methods. The scaling approaches provide qualitatively correct trends of the system behavior over a variety of conditions. The integral equations and the SCF approaches are more ample in the sense that they provide information (e.g., details on the structure of interfaces) which scaling approaches do not. Perhaps, the most accurate are the DFTs that besides an insight into the structure also provide an elegant way to study phase transitions in the systems. Finally, one can classify the models into the ones with or without electrostatic interactions. The systems involving tethered molecules have been also extensively studied using computer simulation methods.

Scaling approaches. A few comments
A system to which the below described theoretical methods apply, is a set of long flexible polymer chains which are anchored on a grafting substrate through a special single end group. In the majority of cases in question, the substrate is taken as flat. The degree of polymerization determines the "chain length" M . On the other hand, the important parameter of the system is the grafting density ρ C , i.e., the fraction of grafted sites per unit area. Usually, this system is put in contact with a solvent. These three parameters (M , ρ C and solvent quality) determine several regimes of behavior of polymer layer having distinct static (structural and thermodynamic) and dynamic properties. The systems of this type are intrinsically characterized by a set of mesoscopic length scales which makes the polymer layers under different conditions amenable to the scaling analysis. Of course, scaling relation imply that all characteristic lenghts are much larger than any atomistic scales, e.g., the range over which the grafting substrate induces density oscillations or the persistence length.
It is commonly assumed that all other segments of chain molecules, besides grafted end group, are not attracted to the surface. Within the scaling theories, the structure of a grafted polymer layer is characterized by the average height and the values of gyration radius along selected axis of grafted molecules. If the grafting density is very low, different coils do not overlap and practically do not interact. Each chain occupies roughly a half-sphere with a radius comparable to the Flory radius for a polymer coil in a good solvent, such that values of gyration radius along selected axis in the z-direction normal to the surface R gz and in the directions perpendicular to the z-axis R g⊥ differ only by their prefactor in the power law relating size to the chain length [1,2], with ν ≈ 3/5, M is the number of segments forming a chain. This behavior is the so-called "mushroom regime".
For denser grafted layers, the chains overlap close to the grafting surface, and are bound to stretch away from the surface, like the bristles in a brush. This situation is quite different from a typical behavior of polymer chains in a solution. If the grafting density ρ C is high (the polymer brush regime) several characteristic length scales appear. According to the scaling approach by Alexander [3] each grafted polymer chain is considered as a linear array of M "blobs", M = M /n blob , where n blob is the number of chain segments per blob (only segements of a single chain are inside a blob), n blob ∝ ξ 1/ν g (ξ g = ρ −1/2 C is the average distance between the grafting points or the blob diameter).
The theory by Alexander, see reference [3] for details, yields the following scaling law for the tethered polymer layer height, h, h ∝ ρ (1/ν−1)/2 This is a very important result showing that the brush height at a chosen fixed grafted density, is proportional to M , while the unstretched chain dimension, R g , grows slower. Consequently, the properties of chains under such conditions may be expected to differ much from unstretched chains in a solution. This scaling argument implies a step-like density profile of the segments. It is difficult to expect a such situation for more realistic models. In fact, entropy favors a nonuniform blob picture. A nonuniform stretching of chains can occur, some chains can turn backwards and some segments can be located anywhere within the brush (of course, there is a maximum distance from the surface within which the segments can be found). Due to nonuniformity, with an increasing distance from the surface, ξ g grows from the value ρ −1/2 C assumed by Alexander, to larger values.
The scaling approaches do not take into account the interactions in the system explicitly. Consequently, specific correlations between particular species of importance in various phenomena cannot be captured in such framework.

Integral equations
Development of theoretical approaches in which intermolecular potentials are the input to the theory is particularly important for discussing phase transitions. One of the first theoretical approaches to describe uniform, as well as nonuniform polymeric systems that explicitly takes into account interaction potentials between all the species was based on the integral equation theory of fluids. Common integral equation theories of atomic fluids were extended to describe polymers by introducing the so-called Polymer Reference Interaction Site Model (PRISM) [5][6][7] that generalized the Reference Interaction Site Model (RISM) by Chandler [8]. For the case of a fluid of homopolymers built of M indentical sites (segments), the PRISM integral equation in the Fourier space reads, where ρ M is the site density (ρ M = M ρ, with ρ being the density of molecules),h(k) is the Fourier transform of the site-site total correlation function,c(k) is the Fourier transform of the site-site direct correlation function. Since all the sites along a chain molecule are equivalent, it follows thath(k) =h α,γ (k) and is the single chain structure factor consisting of Fourier transforms of the site-site intramolecular correlation functions. ω α,γ (k). A generalization of equations (3)-(4) to the cases of mixtures of molecules with different types of sites is available in the literature [5]. After introducing closure approximations, the system of the PRISM integral equations can be solved numerically, yielding the correlation functions that contain a complete information about the microscopic structure of the molecular system. By computing the integrals of the correlation functions, one can obtain thermodynamic properties. The PRISM-based approaches have been used to study a type of model systems involving small solid particles of different geometry (e.g., spherical), coated with polymer chains, see e.g., closely related works [7,[9][10][11]. One recent interesting example is due to Jayaraman and Schweizer [12][13][14] who used the PRISM theory to study model polymer nanocomposites consisting of nanoparticles (called "fillers") with a few grafted chains embedded in a polymer "matrix" (solvent or melt). Coated nanoparticles with a relatively small and controlled number of grafted polymers in this type of mixtures are of interest for various applications. Both the microscopic structure and the phase behavior of such complex systems are sensitive to the effective interactions between species. Those can be tuned by chemical means: synthetic routes of laboratory preparation determine the attraction strength between bare nanoparticles, tether grafting density and matrix chain length. Other important factors are the total packing density and mixture composition, along with the diameter ratio between nanoparticles and polymerizing monomers belonging to grafted chains and to the polymer solvent. We would just like to mention three scenarios leading to the development of different structures.
In general terms, there is a competition between the assumed nanoparticle-nanoparticle attraction, steric repulsion, i.e., excluded volume effects between grafted tethers, and filler-homopolymer correlations. The matrix subsystem can induce a depletion attraction between fillers. Under certain circumstances, the tether-tether repulsion can sterically "screen" the attraction between nanoparticles and if the depletion type entropic attraction is weak, then the fillers form a well-dispersed fluid of a filler in the homopolymer matrix. If the intertether repulsion is weak and the attraction between fillers is strong enough, then the macroscopic phase separation can be induced by depletion matrix effects in addition to bare attraction between nanoparticles. In this case, the matrix-rich and the tethered filler-rich coexisting phases can be formed. Still another, "intermediate" physical situation occurs if nanoparticle-nanoparticle interactions are rather strong and a sufficient steric stabilization due to tether repulsion exists to preclude macrophase demixing, then aggregate formation of microphase-type ordering can take place. So far the PRISM theories have not been applied to describe situations where the grafting occurs on planar surfaces, or, in general, on the surfaces with linear dimension much larger than the size of a polymer segment.

Self-consistent field approaches
One of the popular approaches for the systems of our interest is the self-consistent field theory (SCFT), developed by Edwards [15], DiMarzio [16], Vrij [17] and Helfand and Tagami [18,19]. The key idea behind SCFT is to approximate an ensemble of interacting polymers by the system of mutually non-interacting polymer molecules in an effective, position dependent field. The SCFT methods have been successful in predicting adsorption of polymers, tethered polymers, polymer blends and copolymers, but they are intrinsically incapable of capturing fine details of the system structure at short-range scale [20][21][22]. Due to space limitations, we would like to restrict the discussion of the method to the statement of the problem rather than to the enumeration of several principal results.
For illustrative purposes, we begin with a general statistical mechanical arguments and consider the Gaussian chain model. The system is a blend of N λ kinds of linear homopolymers. The length of the homopolymers of a given kind λ is M λ . For Gaussian chains, the potential energy of the system is [23] βV pot = λα the sum over αλi is over all the species, over all the particles of the given species and over all the sites of a molecule of a given kind and l λ is the corresponding Kuhn length; β = 1/kT , u λλ and v λ are the segment-segment interaction potential and the external field, respectively. The grand partition function for the model in question is, where µ λ is the chemical potential per chain of species λ, and Λ λ = h/(2πm λ kT ) 1/2 is the thermal de Broglie length for monomers of species λ. It is convenient to define the effective field function W λ (r) via the relation where the density operator for the species λ has been introduced Then, the partition function (6) can be rewritten as The segment density profile is obtained from Ξ by using the definition ρ λ (r) = 〈ρ λ (r)〉 = δ ln Ξ δW λ (r) .
Then, the free energy follows from the Legendre transformation [23]. In general terms, the density-field relation (10) can be inverted analytically for an ideal system to yield δW λ (r) as a function of segment density profiles ρ λ (r). Such a task can be realized for the Gaussian chain ideal model that permits to reach the principal equations of the self-consistent field method. Tang and Freed [24] have chosen the model with a vanishing pair interaction, u λλ = 0, as the ideal system such that the potential energy given by equation (5) contains solely the connectivity of chains contribution as well as the term coming from an arbitrary external field. The partition function of the model factorizes as follows: where Q λ denotes the partition function for a single Gaussian chain belonging to the species λ under the effect of the arbitrary external field W λ . The density profile of segments can be easily obtained, such that

12601-5
It is more convenient to develop the following derivation presented below using a continuous chain notation, or, in other words, substituting sums over monomer indices by integrals. The single chain partition function, Q (the subscript λ is dropped to simplify the notation) reads where G(r, r ; M , 0|[W ]) is the unnormalized end-vector distribution function, given as The chain contour variable τ changes between 0 and M . In general, the end-vector distribution function G(r, r ; τ, τ |[W ]) gives the number of chain configurations that start at r , τ and end at r, τ. The path integral representation of the end-vector distribution function given by equation (15) may be alternatively transformed into the partial differential equation that describes a diffusion in an external field W (r). This type of equations serves as the starting point in several research works performed by using the SCFT, see e.g., [9,25]. Additional boundary conditions can be used for a set of specific problems that include confinement effects. It is worth mentioning that the inversion of a field-density relation can be explored in a wider context for systems of interacting chain molecules. The present simple illustration intrinsically leads us to the explanation of the density functional approaches.

Density functional approaches. Generalities
The powerful group of methods to describe nonuniform complex fluids is built based on the ideas of DFT for simple fluids [26]. Perhaps, the first formulation of the DFT for polyatomic fluids was given by Chandler et al. [27,28]. They extended the DFT of nonuniform simple fluids to the systems composed of polyatomic species, considered in the framework of the RISM approach [29,30]. By using Legendre transforms, these authors proved the existence of a free energy functional with the local densities referring to the interaction sites rather than to molecular coordinates. The methodology proposed by Chandler et al. [27,28] retains mathematical simplicity of the traditional theory for atomic fluids, but it requires the "site-site" direct correlation function as an input.
Another class of DFT approaches for polymeric fluids was proposed by Woodward [31], who used the concept of a weighted density approximation, common in the DFTs for simple fluids [32]. Woodward applied the theory to a fluid of particles made of tangentially bonded hard spheres and borrowed the construction of weighted densities from DFT for hard spheres. He obtained the expression for a nonlocal free energy functional that depends on the segment density only. This approximation was also used to investigate a mixture of polymers and hard spheres in slit-like pores. The approach by Woodward was extended to polydisperse polymeric systems [33] with the Schulz-Flory mass distribution, and to infinitely long polymers [34]. At present it is clear that that the weighing procedure of this method is crude and requires improvements. However, a weakness of the theory is due to the presence of an adjustable parameter to account for complex correlation effects.
An important group of DF approaches has been constructed by using the thermodynamic perturbation theory (TPT) by Wertheim [35][36][37]. The earliest theory of this type was proposed by Kierlik and Rosinberg [38][39][40][41]. This approach applies the intramolecular distribution function for pairs of bonded atoms to evaluate the functional of the free energy of the system and a more sophisticated procedure for obtaining the weighted densities.
The theory by Kierlik and Rosinberg was refined by Yu and Wu [42], who introduced a new concept for the treatment of intramolecular degrees of freedom. The novelty is in the use of the fundamental measure theory to accurately evaluate the free energy of a hard sphere reference system and pair cavity 12601-6 distribution function at contact necessary in the TPT formulation. The approach developed by Yu and Wu was used in several theoretical studies of inhomogeneous polymeric systems including surface phase transitions [43], phase equilibria in confined geometry [44][45][46] and effective interactions [47]. The Yu and Wu functional permitted to derive an analytical expression for the polymer-wall interfacial tension [48], to study star polymers [49][50][51], rod-polymer mixtures [52], and hyperbranched polymers [53]. The approach by Yu and Wu has been widely used to study the systems involving tethered chains [54][55][56][57][58][59][60][61][62][63][64][65][66].
Making a short overview of the DF approaches used in the description of nonuniform fluids involving chain molecules, one should mention a quite accurate TPT-based DFT approach by Chapman et al. [67][68][69]. It concerns the associating fluids that can form different complexes and chains, in particular. The DF formulation for associating monomers is supplemented by appropriate mass action laws. The limit of complete association is imposed at the final stage of the procedure, in contrast to the method by Yu and Wu. On the other hand, associative approach is useful to describe the bonding of chain molecules with a solid surface.
It should be noted that there exists a link between DF methods and the SCFT for systems of polyatomic molecules. Several DFTs can be recast to the SCFT form [23,24] with a different degree of difficulty. Generally, the so-called DFT-SCFT approaches rely on a combination of computer simulations of a single chain in an effective position dependent field that results from nonlocal DFT [70][71][72][73][74][75]. The advantage of the DFT-SCFT approaches is that practically all global and local properties of polymers at interfaces, such as the center-of-mass profile or the end-to-end vector profile, can be calculated and analyzed [76][77][78][79]. In particular, those methods can be used to take into account the shrinking or the swelling polymers at the surface. So far, such calculations have been carried out only for bulk systems using PRISM integral equation approaches [80][81][82].
Quite recently, Bryk and MacDowell proposed a new hybrid self-consistent field DF approach to describe the solvation effects for polymers at interfaces [83]. Using DFT of Kierlik and Rosinberg [38][39][40][41], they showed that when the second-order terms of Wertheim thermodynamic perturbation approach (TPT2) [35][36][37] are taken into account, the DFT-TPT2 is capable of capturing the effects of solvation of tethered chains. Theoretical predictions were tested with respect to Monte Carlo simulations for moderate chain lengths [83]. The results for the end-to-end distances of chain molecules in the solvent in the bulk phase and at the solid surface were in a reasonable agreement with simulations.

Density functional theory
As we have already mentioned in the introductory section, different formulations of DFT are available in the theory of inhomogeneous complex and polymeric fluids. In this text, however, we restrict ourselves to solely one, frequently used version, for illustrative purposes. It contains elements of other approaches and in the course of presentation we will address the issues related to other DFTs.
Let us consider the fluid (indexed as 1) in contact with the surface covered with the end-grafted chains, P , [54][55][56][57][58][59][60][61][62][63] (an extension to fluid mixtures is straightforward). The chain consists of M tangentially bonded segments. Although other models involving different segments can be considered, we restrict our attention to the case where all the segments are identical.
The connectivity of segments within a chain is provided by the bonding potential, V B , defined as follows: where R ≡ (r 1 , r 2 , · · · , r M ) is the vector specifying the positions of segments and σ (P ) is the diameter of the segments. It is assumed that all the segments and fluid molecules interact with each other and with solvent molecules by means of (12-6) Lennard-Jones potential, where ε (kl ) is the strength of interactions between the species k and l and σ (kl ) = σ (k) + σ (l ) /2, (k, l = P, 1).

12601-7
The chains are pinned to the surface. The interaction of their first segments with the surface is given by where z is the distance from the surface and C is a constant. The remaining segments of the chains, j = 2, 3, . . . , M , as well as the fluid molecules interact with the surface via a hard-wall potential, or, depending on the model used, via the Lennard-Jones (9-3) potential, In the above equation, ε (k) s is the strength of interaction and z 0 = σ (k) /2. The geometry of the system can be changed. One can study a slit-like pore, then the potentials (19), (20) and (21) must be modified, cf. [57,58].
In order to present the theory, let us introduce the notation ρ (P ) (R) and ρ (1) (r) for the density profiles of chains and fluid, respectively, and define the segment densities ρ (P ) s j and the total segment density of chains, ρ (P ) s , via the relation, To simplify the equations, we use the notation ρ (1) s (r) ≡ ρ (1) s1 (r) ≡ ρ (1) (r). The system is considered in the grand canonical ensemble. The thermodynamic potential, Ω, is as follows: where F ρ (P ) (R), ρ (1) (r) is the Helmholtz free energy functional and µ is the chemical potential of species 1. The functional Ω is evaluated under the constraint where A is the surface area and ρ C is the surface density of grafted chains. The free energy functional is divided into ideal, F id , and the excess, F ex , terms. The excess terms are due to the hard-sphere interactions (hs), due to the connectivity of segments of chain molecules (c) and due to attractive interactions between all the species (att), The configurational ideal part of the free energy functional is as follows: The excess free energy, due to hard-sphere interactions, F hs = Φ hs (r)d r, results from the "white bear" version of the fundamental measure theory. Then, we have [42] In the above, ξ(r) = |n 2 (r)|/n 2 (r). The definitions of scalar, n α , for α = 0, 1, 2, 3 and vector, n α , for α = 1, 2 weighted densities are given in reference [42].
Finally, the attractive contribution is given by the mean-field approximation where u (kl ) att (|r − r |) is defined according to Weeks-Chandler-Andersen scheme, and the diameter is assumed d (i ) = σ (i ) , (i = P, 1).
The equilibrium density profiles are obtained by minimizing Ω under the constraint (24), The final equations for the segment density profile and for the fluid density profile can be found in references [42,57,58]. The bulk fluid in equilibrium with the system of interest is the fluid of Lennard-Jones particles of type 1. The configurational chemical potential µ in terms of the bulk density ρ b of species 1 is where µ (hs) is the excess chemical potential of hard spheres.
The knowledge of the density profiles is crucial for the evaluation of thermodynamic properties. The density profiles provide the excess adsorption isotherms of fluid species, Γ ex , where the limits of integration are determined by the geometry of the system. The knowledge of the density profiles yields the value of Ω and thus permits for the evaluation of the solvation force The solvation force between the surfaces is a measurable property. It comes out from the experiments performed by using a surface force apparatus, see e.g., [84]. Our final comment concerns another important property. Upon fluid adsorption, the average height of tethered chains changes in a nontrivial way. In order to describe this property, we use the normalized first moment, Depending on the physical argument, the coefficient a is taken equal to 8/3 (cf. [85,86]) or a = 2 [87].
The theory presented above permits for several extensions. In particular, it is possible to consider either tethered copolymers [88] or chain layers of different species. Flexibility of bonding between segments can be modified as well, assuming, for example, preferential angles for bonding between certain groups of segments [89]. It is also possible to consider the branched polymers and to change the identity of segments (and their sequence).

Applications of the density functional approaches
In this section, we would like to focus on a few results of DFTs. In general, one needs to be convinced that the obtained density profiles and the values of Ω are accurate for different values of parameters describing the system. The profiles yield the excess adsorption and the values of Ω allow for the investigation of the phase behavior.
We proceed to a discussion of density profiles. The configuration of end-grafted polymers in a good solvent has been well documented. Specifically, the SCF theory by Milner et al. [90] based on the ideas of Alexander [3] and de Gennes [1], lead to conclusions that at a moderate surface grafting density, the chains exhibit parabolic density profiles. This result has been supported by simulation data. We focus on exploring the possibilities of DFTs and the DFT by Chapman et al. [69], denominated as iSAFT (interfacial statistical associating fluid theory). As we have already mentioned, this approach slightly differs from the theory described above by the treatment of bonding of chain molecules to the solid surface.
The first illustration [69] shows the structure of athermal chains with M = 100 hard-sphere segments tethered to a hard wall. The density profile is parabolic (left hand panel of figure 1). At a high grafting density, the deviations from a parabolic shape are not big. In general, this DFT overestimates the density profile values close to the wall and underestimates the profile further from the wall. However, the difference between simulation data [91] and theory is not large. In addition, Chapman et al. [69] showed that the DFT yields predictions similar to the SCF approach [92]. The right hand panel of figure 1 shows the density profiles of chains tethered to a planar surface in the presence of free polymer solution of short chains resulting from iSAFT theory and compares them with simulation results and with the profiles for an implicit solvent. The density profiles for the models with explicit solvent significantly differ from the implicit solvent case. The segments of solvent chains interact via a purely repulsive potential. Consequently, the observed effect is entropic. In an explicit solvent, the solvent molecules compress the chains and the layer height is reduced compared to an implicit solvent N c kind of linear homopolymers case. The total solvent density is high, see the figure caption. The compression increases with an increase of the number of segments of the solvent, showing that penetration of longer (M = 10) free chains into the tethered layer decreases, compared to shorter chains. Thus, the chain solvent effect is to reduce the tethered layer height and to make it denser in the vicinity of the solid surface. The observations from the DFT are in agreement with computer simulations [91].

Phase transitions
After evaluating the density profiles, one can try to describe the thermodynamics of adsorption and to analyze the phase behavior. Such problems were studied by Borówko et al. [62], who considered a monoatomic fluid in contact with modified and non-modified solid surfaces. The fluid-wall interaction was given by the Lennard-Jones (9,3) potential. The surface was modified by short (5 segments) chain molecules. The model assumed the same interactions ε = ε (kl ) between all species, the diameters of segments and of fluid particles were the same and ε s /ε = 8; ε s = ε (P ) s = ε (1) s . The DF approach described in section 3 was used.  A set of adsorption isotherms, at different reduced temperatures, T * = kT /ε, is shown in figure 2 (a). At high temperatures, the adsorption isotherms, Γ * = Γσ 2 , monotonously increase with bulk density, ρ * b = ρ b σ 3 . However, at lower temperatures, they exhibit hysteresis loops. The first-order phase transition point at each temperature was found from the crossing of two branches of the thermodynamic potential, In order to characterize these phenomena it is necessary to construct the phase diagram. Numerical calculations revealed that for the system without tethered chains, the wetting temperature is T * w = kT w /ε = 1.03, whereas the surface critical temperature is T * sc = kT sc /ε = 1.125. Upon introducing even a small amount of chains, one observes a shift of both characteristic temperatures to lower values. At the grafting density ρ * C = ρ c σ 2 = 0.01, the values of T * w and T * sc are T * w = 0.83, and T * sc = 0.92, respectively, see Similarly to the changes of the wetting properties, one can investigate layering transitions. The layering transitions are usually discussed for fluids in contact with strongly adsorbing surfaces. The aim is

12601-11
to explore how tethered chains change the scenarios for layering transitions. In particular, whether it is possible to suppress the layering of a fluid or to eliminate some transitions. However, for some systems, one can also observe new layering transitions as a consequence of chemical modification of a weakly attractive surface by strongly attractive chains. The important factors effecting the changes of the phase behavior are the differences in the interaction energies between the bare surface and the segments of chains and the entropic effects due to the changes in the structure of chains. The parameters such as the number of segments in the chains and grafting density are of importance as well. These factors permit to tune the phase behavior in terms of the formation or destruction of consecutive layers. In the case of adsorption in pores, the modification of the pore walls effects the capillary condensation. In this case, the ratio between the chain length and the pore width becomes a new important parameter. The chains can also enhance or inhibit layering transitions that may preceed the capillary condensation. All these problems were discussed in references [58,59].
If the binding of the terminal segments of chains to the walls of a slit-like pore is not strong, a new kind of phase transition can take place. This transition is connected with the change of the symmetry of the distribution of chains and fluid molecules inside the pore. Upon the adsorption of a fluid, a part of weakly tethered chains can be "wiped out" from one wall and bonded to the opposite one. Consequently, the distribution of the chains inside the pore becomes non-symmetric. As a result, the symmetry of the distribution of fluid molecules with respect to the pore center is also broken. A further increase of the fluid density may induce a reentrant behavior manifested in the recovery of the symmetry. Such changes of symmetry may occur as first-or second-order transitions [64].
The symmetry breaking transitions are also observed for slit-like pores modified with chains whose two ends are permanently attached to the opposing walls, provided that the length of chains is sufficiently larger than the pore width [65,66]. Such chains are often called pillars. In this case, the inner segments of the chain molecules can move towards one of the pore walls. This problem was discussed in references [59,62].
The phase transitions, associated with the changes of the symmetry of local densities can be observed if the chains are much longer than the pore width. For chains of the length comparable to the pore width, such transitions do not occur [65]. Moreover, the symmetry changes are possible if the surface density of the chains is moderate. In the limiting case of a very small grafting density, the phase transformations in modified and non-modified pores are similar. In another limiting case of a very high grafting density, all phase transitions within the confined fluid disappear. The values of the fluid-wall energy parameters are also important. In the case of strongly attractive adsorbing potential, the symmetry breaking is linked with the layering (that takes place at one wall). For weaker adsorbing potentials, the growth of a non-symmetric film at one wall resembles the growth of a film after the prewetting jump. Finally, for a very weak adsorbing potential, the symmetry breaking does not occur. In the case of a strongly adsorbing surface, the first-order symmetry change transitions end up at the tricritical points from which the second-order transition lines depart. There exists a triple point between the "return-to-symmetry" transition and the capillary condensation envelope. In the case of weaker adsorbing potentials, however, the second-order symmetry change line joins the capillary condensation part of the phase diagram at the critical end point. In such cases, a symmetry recovery at lower temperatures takes place during capillary condensation.
The phase behavior becomes even more complicated if the adsorption of a binary mixture inside the modified pore is considered. If the mixture exhibits a demixing in the bulk phase, then demixing transitions inside the pore, or even within particular layers formed upon the layering transitions, compete with the symmetry breaking. This competition may lead to complex phase diagrams that include critical points, lower and upper critical end points, tricritical points, λ-lines connected with demixing as well as λ-lines connected with symmetry breaking and symmetry recovery transitions. Some of these problems were discussed in reference [66].

Solvation force
There are several works focused on the study of the effective interaction between surfaces coated with chains. In particular, Wu et al. [54] studied the solvation force between two surfaces covered with telechelic chains. Telechelic chains are chain molecules with certain functional groups at both ends.
Thus, the end segments of a polymer can bind to a single surface simultaneously to form a loop, or to the opposite surfaces to form a bridge-like structure. The backbone segments in between the two terminal groups are chemically inert. Experimental and theoretical studies show that the interaction between highly stretched telechelic chains is mostly repulsive, in close similarity to the behavior observed for the system of two layers of singly-tethered chains. If the two layers of telechelic chains are put almost at contact, i.e., if their separation distance is close to the double height of an isolated chain, a weak attractive force is observed. The origin and the magnitude of this attraction has been discussed by different authors [93][94][95]. Similarly to the DFT by Chapman et al. [67][68][69], the approach by Wu et al. [54] involves a parameter describing the strength of an attractive interaction between terminal segments and the surface, surface-adhesive energy, ε Aw . In figure 3, taken from [54], the solvation force between two identical surfaces covered by telechelic chains is shown. The dependence of the solvation force on the value of ε Aw and on the number of segments, M , was studied. The segment packing fraction in the bulk phase of chains was rather low, ρ * b = 0.05. If the surface-adhesive energy increases from ε Aw /kT = 1, to 8, 12 and 15, the solvation force exhibits changes from a short-range attraction to a long-range repulsion. At the lowest energy ε Aw /kT , the solvation force is similar to that for nonadsorbing polymers, ε Aw /kT = 0; the attractive short-range part of the force is due to the entropic depletion of chains from the surface. At the highest value of ε Aw /kT , the repulsive force arises from the interaction between tethered telechelic chains. The strength of repulsion depends on the energy ε Aw /kT . The magnified view of the solvation force shown in the inset indicates that in a nonadsorbing regime (small ε Aw /kT ), the theory leads to entropic depletion, as well as to a weak, secondary repulsion between surfaces. If the polymers are strongly adsorbed, the theory predicts the existence of a weak attraction between telechelic chains at contact. Moreover, the location of the attraction shows that the chains are highly stretched.

12601-13
DF observations are in agreement with experimental findings and other theories. If the chain length increases, similar conclusions concerning the behavior of the solvation force preserve. A more detailed discussion can be found in reference [54]. To summarize our illustrations in this part, it is worth mentioning that similar successful developments have been used for the description of inhomogenous fluid mixtures, tethered chains of copolymer type and other related models with nonelectrostatic interparticle interactions. We proceed now to another class of systems, namely to the models with long-range Coulomb interaction.

Electrostatic density functionals
An ample set of problems of interest for laboratory and for applications involves fluids of charged particles, as well as tethered polyelectrolytes and polyampholytes. In the case of tethered polyelectrolyte chains, the presence of long-range electrostatic interactions makes the properties of such systems much more complicated than those corresponding to neutral grafted chains. Since the pioneering works by Pincus [96], Balastre et al. [97] and Borisov et al. [98], numerous theoretical and experimental studies have been reported concerning the structure and surface properties of charged grafted chains. Similarly to the case of uncharged chains, the theoretical investigations were primarily based on the scaling analyses and the self-consistent mean-field theories that all rely on a number of simplifications [96,[98][99][100][101][102]. If the number and positions of charges at each chain are fixed, theoretical and experimental investigations elucidate the existence of different brush regimes that depend on the grafting density, degree of dissociation, and ionic strength. However, more recent experiments by Wu et al. [103] and Currie et al. [104] have shown that the predicted scaling relations do not necessarily hold throughout the entire range of the grafting densities and ionic strengths. Recently, Witte et al. [105] developed a new version of a SCF approach from the grand canonical partition function. The theory properly accounts for the local nature of the charge equilibrium and captures the basic behaviors of grafted polyelectrolytes. In agreement with recent experiments, the theory predicts that the scaling of grafted chains height with grafting density can be qualitatively different at intermediate chain lengths from that predicted by scaling arguments. This difference was attributed to the relative strength of electrostatic interactions compared to finite segment size packing constraints. On the other hand, the trend of decreasing chains height with an increasing grafting density predicted by the scaling analysis is recovered for large molecular weight polymers immersed in a solution of a very weak ionic strength.
One should also mention here the SCF approach [106] aimed at predicting the structural properties of polyelectrolytes grafted via one of their ends to solid surfaces which explicitly incorporates the acid-base equilibrium responsible for the charge regulation of the polyelectrolyte groups, as well as the conformations, size, shape, and charge distributions of all the molecular species present. The rest of that theory was compared with experimental observations for the height of the layer as a function of ionic strength for different chain molecular weight and grafting densities, and a good agreement was found.
Similarly to the case of uncharged chains, DFT provides an alternative method to traditional scaling and mean-field theories by taking into account both short-and long-range correlations in a self-consistent manner [107,108]. In order to illustrate the current developments of the DFT for this type of systems, we would like to briefly consider the adsorption of an electrolyte solution on the solid surface modified by grafted chain molecules containing charged segments. Commonly, an electrolyte solution is modelled in the framework of the primitive model, i.e., the fluid is considered as a mixture of hard-spheres interacting via Coulomb potential. The same form of the potential is assumed for the interaction between the fluid species and charged segments of the tethered chains, where the superscripts η, α = P, 1, 2, and Z α i is the valency of ions or of the charged segments. The solvent is treated as a continuum with a given relative dielectric permittivity, ( 0 is the permittivity of the vacuum).
The interaction between the ions and the wall, as well as between the charged segments of the chains and the wall, contains the electrostatic term, in addition to non-electrostatic interaction. Here, Qe is the surface charge density on the wall and l B = e 2 /(4πkT 0 ) is the Bjerrum length. The bulk densities of ionic species of the solution must satisfy the electroneutrality condition Z (1) The electroneutrality holds for the entire system and the charge of tethered chains is neutralized by charges of counterions.
The distribution of all charges in the system satisfies the Poisson equation, where Ψ(z) is the electrostatic potential. The solution of this equation requires boundary conditions. For a fluid in contact with a single (modified) wall at z = 0, we can specify the value of the electrostatic potential at the wall, setting V 0 = Ψ(z = 0). The wall charge follows then from the electroneutrality condition, 1 ρ (η) (z). One can work either under the condition of a constant Q or V 0 .
Similarly to the previous section, the principal issue is in the determination of the density profiles of fluid species and of tethered chains segments by minimization of the grand thermodynamic potential. In the present problem, the excess grand potential is given by In contrast to section 3 the free energy F also contains the term due to the coupling between electrostatic and hard sphere interactions [109,110]. It is given by The summation in equation (41)  are short-range parts of electrostatic two-particle direct correlation functions, that can be taken from the analytic solution of the mean spherical approximation for a restricted primitive model of electrolyte solutions and weighted correlation approach by Wang et al. [109,110]. Another difference between the systems involving charged and uncharged segments is that the free energy term describing the connectivity of chains is modified. According to Wertheim's theory [35][36][37] for chains built of tangentially jointed hard spheres, the free energy contribution F c (cf. section. 3) involves the contact value of the cavity function of the reference system of hard spheres, y (P ) hs . For uncharged systems, a key assumption of the DFT is that the contact value of the segment-segment cavity function is represented by that corresponding to a fluid composed of hard spheres, y (P ) hs . In the case of chains built of charged segments, an approximation for the cavity function includes the effects of screened segmentsegment electrostatic interactions and takes the form [111,112], In the above, Γ = 1 + 2κσ − 1 /2σ, T * = σ/l B is the "electrostatic" reduced temperature and

12601-15
is the square of the "local" screening parameter [109,110]. The application of the above approximation for chains with different charge sequences under different thermodynamic conditions leads to qualitatively correct results [112].

Selected applications of the electrostatic density functional approach
In general, tethered polyelectrolytes respond to external stimuli such as the properties of electrolyte solutions, density, ionic strength, temperature, see, references [113][114][115][116], as well external electric field. There exist a vast literature dedicated to experimental, simulation and theoretical studies of such responsive chains. Among various external stimuli, electric field is the one of much interest because it can induce changes of the tethered layer height [114,117,118] and consequently yields changes of, e.g., permeability, transport properties and electric conductance of nanochannels [119].
Upon putting an ionic fluid in contact with the modified surface, the height of the tethered polyelectrolyte changes depending on different factors. A set of results concerning these trends from the DF approach by Jiang et al. [55] is given in figure 4. In this particular case, the polyelectrolyte grafted at the uncharged solid surface consists of univalent negatively charged spherical segments, whereas the fluid (salt) is a mixture of cations and anions of the same valency. For simplicity, the diameters of ions and segments are assumed to be equal. At a given chain length, the height of polyelectrolyte layer practically does not depend on the salt density in the interval of low salt densities. This is the so-called osmotic-brush regime and this type of behavior was first established theoretically, see e.g., [96] and confirmed experimentally [97]. By contrast, at higher salt concentrations, the results shown in the figure are in accordance with the salted-brush regime, in which the brush thickness obeys the power law, established previously in references [96,98]. In the unifying picture for a wide range of salt concentrations, one can say that polyelectrolyte chains undergo a smooth transition from "osmotic sticks" to salted coils. The DFT reproduces the previously established scaling laws for long chains. Moreover, for shorter chains, e.g., for M = 20, the DFT provides very good estimates for the tethered layer height on salt density in close similarity to other approaches.
Another result concerns the dependence of the height of chains on the grafting density. Although the grafting density in the previous example was low, this value still yields the brush rather than mushroomtype structure, as it can be deduced from the values in the figure. The grafting density parameter effects the height of the chains in a nontrivial manner. Namely, in the osmotic regime, the height of the chains increases slightly with an increase of grafting density showing deviation from the scaling law behavior. Such a growth is in accordance with computer simulations and recent experimental findings [120,121].
This "nonlinear osmotic" regime (i.e., when the concentration of counterions inside the tethered chains is lower than in the bulk fluid) is well described by DFT. Moreover, the slope of the curve changes with the grafting density, or, in other words, there is no simple scaling relation between two variables of the figure at low salt concentration. In the "salted-brush" regime, the height of the chains grows slowly until the separation between the grafting points of polyions becomes small, such that at a higher grafting density, the height of the chains growth obeys a quasi-power law. In fact, the emerging picture is richer and gives the details that are not intrinsically provided by simple scaling laws.
Previous examples have concerned the polyelectrolytes tethered at uncharged surfaces. However, there is an increasing interest in developing nanofluidic systems that can control various properties, including transport of ions and molecules in aqueous solutions in charged nanochannels. The next example shows the effects of non-zero electrostatic potential at the surface on the distribution of chain segments and of ions in the interfacial region. One would expect that the electric double layer will be substantially affected by the presence of tethered chains. Our presentation concerns the primitive model of electrolyte. Moreover, we assume that tethered chains consist of uncharged segments of the same diameter as fluid ions. Two walls covered with tethered chains constitute a slit-like pore. We apply the theory briefly described in the first part of this section. The walls are charged, and the charge is determined by the electrostatic potential on the wall V * = eV 0 /kT . One of the properties of interest is the differential capacitance. It is defined as the derivative of the surface charge with respect to the applied electrostatic potential, V * , Recently there has been much interest in the dependence of a differential capacitance on the pore width [122][123][124][125][126]. This interest has been inspired by a possibility of developing new electric double layer capacitors. In particular, it is quite challenging for the theory to provide interpretation of experimental observations concerning an anomalous growth of capacitance with a decreasing pore width for very narrow pores. A few examples of the dependence of the reduced differential capacitance on pore width are given in figure 5. Reprinted with permission from reference [128]. Copyright AIP.

12601-17
In the absence of tethered chains, the dependence of C * d on the reduced pore width, H * = H /σ, exhibits oscillations, with two well pronounced maxima at H * ≈ 1.3 and 3.3, (V * = 0.5) figure 5 (a). The maxima are separated by a trough. In the presence of tethered chains made of six uncharged segments (this system is abbreviated as N6), the magnitude of oscillations decreases, but at a low grafting density, ρ * C = 0.05, the capacitance curve still preserves the shape similar to the case without chains. At higher grafting density, ρ * C = 0.1, the first maximum is washed out and only the reminiscence of the second maximum preserves.
Even a relatively small amount of chains substantially alters the electric capacitance in narrow pores compared to a pore with bare walls. Since the segments are uncharged, the observed effect of chains is due to excluded volume effects. Figure 5 (b) shows the results for the pores modified by grafted polyampholytes, i.e., copolymers consisting of positively and negatively charged segments (in some models uncharged groups are considered as well, see e.g., [127]). The segments with charges of different sign can be distributed randomly or can form certain patterns. In particular, the charges of the same sign can be arranged in blocks. The results displayed here have been obtained for the chains whose net charge is zero, i.e., the number of positively and negatively charged segments within a chain is the same. Two polyampholyte chains were considered [128], namely built of 6 and of 10 segments, whose 3 or 5 initial segments bear positive unit charges, while the following segments bear negative unit charges. These systems are abbreviated as P3M3 and P5M5, respectively. The majority of the results in figure 5 (b) are for a low grafting density, ρ * C = 0.025. Higher grafting density suppresses oscillations, similarly to the case of uncharged segments. The electrostatic potential at the walls is low because pronounced oscillations appear only at low voltages, as in the case of non-modified pore walls [58]. Dashed lines in figure 5 were obtained at V * = 3. Solid lines, however, were evaluated at V * = V PZC , the potential of zero charge.
Note that for the system involving polyampholytes, the value of the potential of zero charge differs from zero, even for chain particles whose net charge is zero [130]. Of course, the value of V PZC changes with the pore width. In the case of non-modified pores, the first (and the most significant) maximum appears at the distance close to H * ≈ 1.3. The presence of a brush significantly diminishes this maximum. The maximum still preserves for the brush P3M3, but vanishes for the brush of longer chains, P5M5.
In the case of a pore with non-modified walls [129], the the distribution of charges inside the pores, and, consequently, the dependence of Q * on V * , results from the possibility of the development of a number of layers of ions at both walls and from the interference between the layers formed at the opposing walls. For pores having walls covered with brushes whose segments are charged, there exists a constant amount of charges inside the pore, due to charged segments. Their distribution is not "free", but is restricted by the existing bonds between the segments. This "background" charge distribution reduces the freedom of changes of the distribution of ions and, consequently, variations of Q * with V * become smaller. This is reflected by a smaller (in comparison to non-modified pores) amplitude of the function Another interesting question is the effect of the brush on the dependence of the double layer capacitance on the voltage (or the charge) at the wall. It is known that for non-modified surfaces, the shape of this curve changes from a single-hump ("bell-like") to a double hump ("camel like") if the bulk density changes. For non-modified surfaces, the appearance of single, as well as double humps was found in Monte Carlo simulations by Fedorov et al. [131,132], and in DF calculations by Wu et al. [133] and [134]. This behavior was experimentally observed in the systems that involve surfaces free from grafted species [135,136].
The calculations carried out in reference [128] have revealed that for narrow pores modified with grafted polyampholytes, even three humps can appear on the dependence of the capacitance on the electrostatic potential on the wall and that the shape of the capacitance curve depends on the grafting density. In general, in the presence of tethered species in slit-like pores, the shape of C * d (Q * ) changes from a camel-like to a bell-like, if the bulk RPM fluid density changes from low to moderate and high. The value of crossover bulk density depends on the grafting density of the chains, their length and on the charge sequence along the chains. It is possible to alter the shape of C * d (Q * ) curve from double to single hump by changing the grafting density, if the bulk fluid density is adequately chosen. The symmetry of C * d (Q * ) curves with respect to Q * = 0, observed for several model brushes of uncharged segments, is destroyed if one explores the brushes made of the chains having charged segments.
To summarize, these few examples show the capability of the DF approaches to describe the models of complex systems comprising a solid surface with end-tethered polyelectrolytes or polyampholytes and electrolyte solutions. There is much room for improving the modelling and theory. One of the principal problems is an adequate description of the aqueous solvent and its dielectric properties, as well as to reasonably capture the screening of ion-solvent interactions. The inclusion of orientation degrees of freedom of species into the DFTs is, therefore, necessary. In general, we expect that the systems described above will find their place in the research aimed at optimizing the energy density of supercapacitors, see e.g., [137].

Computer simulations
Computer simulation is considered as an intermediate method bridging up the experimental and theoretical studies. Theoretical DF approaches have already been discussed above. Concerning the experimental studies, one should note that recent technological progress in the fabrication of modified surfaces involving tethered chains is not adequately followed by the advances in the methods for their experimental analysis [138]. The existing techniques such as optical, X-ray and neutron scattering, operate on the meso-and micro-length scales and cannot be informative enough for the systems rapidly changing on the nano-scale (as in the case of nano-patterned surfaces). On the other hand, the Atomic Force Microscopy (AFM) can be invasive [139], and can lead to non-negligible changes of the material structure [138]. Other difficulties are connected with a precise control of polymer synthesis (problems with polydispersity, defects, purity, etc.).
The benefit of computer simulations is that they provide access to the positions (and in some methods to the velocities) of individual particles enabling a high resolution insight into the structure and internal dynamics of the particles constituting the material. As a result, the relevant microscopic mechanisms of particular phenomena can be clarified or even predicted in the way inaccessible for theory and the experiment. This level of resolution is achieved at the expense of simplification of the interaction potentials (often being empirical) and of a finite size of the simulated system. Most simulations of the systems involving tethered chains can be split into two categories [140][141][142]: the (semi-)atomistic (the particle represents an atom or a group of atoms involving hydrogens, e.g., CH 2 ) and mesoscopic (the particle represents a fragment of a polymer chain or a collection of solvent molecules). Both approaches can be thought of as mutually complementary. While the (semi-) atomistic modelling allows one to address the real chemistry of the system (and the effects like hydrogen bonding, entanglements etc.), it lacks a possibility to cover large length-and time-scales of the phenomena under study. The mesoscopic modelling is capable of doing that, but at the expense of introducing a certain level of coarse-graining focused on chemical (in)compatibility of fragments and the effects of their microphase separation. For further details of various particle-based simulation techniques for the systems involving tethered chains see the reviews [140][141][142][143][144].
DPD is an off-lattice mesoscopic simulation technique which involves a set of particles moving in continuous space on a discrete time scale. Particles (beads) represent whole molecules, fragments of the polymer chain or fluid regions, rather than single atoms, and atomistic details are not considered. In other words, the internal degrees of freedom of particles are integrated out and instead of real interparticle interactions one introduces effective pairwise forces that depend on bead-bead distances. The total force acting on i -th bead from the interaction with all other j -th particles is given as a sum where F C i j is a conservative force that identifies polymer beads, while the dissipative, F D i j , and random, F R i j , forces form a thermostat that keeps the temperature of the system constant. All the details of the

12601-19
method (e.g., justification of introduction of effective forces and explanation how to derive them) are given elsewhere [150][151][152]. A key property of the forces in DPD simulation is that they conserve the momentum locally, so that hydrodynamic modes of the fluid emerge even for small particle numbers. The parametrization of conservative forces is density and temperature dependent and can be linked to the Flory-Huggins parameter for the case of mixtures (see, e.g., reference [154]). The gain of such a way of coarse-graining is twofold. Firstly, the number of moving entities is essentially reduced. Secondly, the effective forces in DPD are much "softer" than the atom-atom ones. This allows one for using much larger timesteps in the DPD simulation, in comparison with (semi-)atomistic MD discussed below. The conservation of total momentum ensures a correct hydrodynamic behavior of the systems (i.e., the system satisfies the Langevin hydrodynamic equations). As a result, the method permits to access much longer time and length scales and ensures a faster micro-phase separation of non-miscible species, as compared with conventional MD simulations. Simulations of polymeric fluids in volumes up to 100 nm in linear dimension for tens of microseconds are now common [151].
The DPD modelling of confined fluids, or systems with walls, requires a special care. As it was discussed by Visser et al. [155], the central problem is how to construct a solid wall for mesoscopic modelling. At this level of modelling, three boundary conditions must hold for a solid wall: (i) impenetrability (particles are not allowed to enter the wall), (ii) the wall should not artificially affect (i.e., by imposing an artificial structure) the fluid properties, and, in the case of flows, (iii) no-slip; the wall should impose velocities in accordance with Langevin hydrodynamics. Previous studies succeed to implement impenetrability and no-slip boundary conditions. They usually used several layers of frozen particles residing on a surface [156][157][158][159] to ensure impenetrability. This, however, induces the ordering near a surface, which propagates into a substantial portion of the system. This effect has a physical justification for the atom-based simulation when one mimics the crystal-like structure of the wall, but should be interpreted as artificial for the meso-scale type of modelling.
Quite recently, a new approach was introduced for modelling the walls in DPD studies [160]. It was assumed that, if the i -th bead with coordinates r i = {x i , y i , z i } approaches the surface, e.g., the plane z = 0, it interacts with an "imaginary" bead that has the coordinates r j = {x i , y i , 0}. The total wall-bead interaction has all three contributions given by equation (45). Repulsion of the surface can be tuned by adjusting the parameters of conservative force, F C . The condition of total momentum conservation for the wall-bead interaction is applied by the following procedure. The force acquired by each wall due to the interaction with the approaching beads is accumulated during the evaluation of forces. Then, it is evenly distributed over all the beads inside the simulation box. The effective surface, therefore, is structureless and, thus, it does not introduce any artificial perturbation of the fluid structure near the wall. The details concerning DPD simulations can be found in references [151,156,[160][161][162][163][164].
Let us now briefly summarize the most interesting DPD results for the systems involving tethered polymers. One of the first DPD studies of polymer-coated surfaces was carried out by Gibson et al. [165], who investigated the adsorption of colloidal particles on a surface modified with pinned polymers. In particular, it was found that the adsorption of colloidal particles was smaller when the size of the polymer relative to the colloidal particle and the density of the polymers increased. Moreover, the adsorption reduces when the polymer is well solvated.
Malfreyt and Tildesley [163] simulated the brushes tethered at two walls of a slit in a good solvent. The density profiles of polymer segments across the pore were found to exhibit parabolic shape, and diffusion along the pore axis was significantly greater for the solvent particles in the middle of the slab than in the polymer region, since the solvent molecules were trapped within the entangled polymer. Later, Goujon, Malfreyt, and Tildesley [156] extended DPD method to the grand canonical ensemble and studied the compression of grafted polymer brushes in a good solvent. They also studied the compression of polymer brushes under shear and evaluated the friction coefficient as a function of compression, shear rate and properties of the solvent [164]. The constant and oscillatory flow of a fluid between two plates with tethered chains was simulated by Wijmans and Smit [166].
Pastorino et al. [167] investigated a short-chain melt between two brush-covered surfaces. The equilibrium situation, as well as the one under shear flow were explored. The polymers of both brush and melt were assumed identical. They studied the interdigitation of the melt and the brush and also considered the orientation of bond vectors on different length scales, as well as radii of gyration, end-to-end vectors of free and grafted chains, and in the case of flow simulations, the velocity profiles. In the second work, Pastorino et al. [168] used the Gibbs criterion to localize the brush-melt interface and analyzed its equilibrium fluctuations in terms of a capillary wave Hamiltonian augmented by an elastic term that accounts for the deformability of the brush. However, one should stress that the simulation method used by Pastorino et al. [167,168] was different from the standard DPD scheme [152]. Namely, they assumed that the interactions between the polymer segments and polymer segments with surfaces were of the Lennard-Jones type, and they performed molecular dynamics simulations with DPD thermostat.
The effects of chain architecture on the structure of tethered rod-coil polymer layers was studied by Li et al. [169]. When immersed in a selective solvent for the coil blocks, rod blocks tend to form aggregates. Linear and Y-shaped polymers exhibited a similar aggregative behavior, but comb-like brushes were found to possess more diverse aggregative behavior compared to linear brushes. Surface structures with aggregates taking form of cones, cylinders, or layers of spheres were found. The behavior of grafted binary polymer brushes with compatible components in the case of different chain lengths was investigated by Xue et al. [170]. They observed layered structures parallel to the surface indicating "phase separation": short chains were suppressed in the layer adjacent to the surface, whereas longer chains were to much extent stretched. By slightly changing the solvent selectivity in preference of shorter chains, inversion of the layered structure was found.
There have also been performed several DPD studies of the brushes tethered to non-uniform surfaces, see e.g., [138,160,161,171,172]. The systems involving brushes attached to some parts of the surface, the so-called structured brushes, are potentially very important in the manipulation of fluids at very short length scales. A prototype case consists of a solvent in contact with a surface, or confined between two identical plane-parallel substrates, decorated with stripes of tethered chains that alternate periodically in one direction (say, X ) and are infinite in the other transverse direction (Y ). The performed simulations indicate that heterogeneity of tethered layers has a great impact on the structure of confined fluid and on thermodynamic and dynamic properties of the systems. In particular, Ilnytskyi et al. [160] employed DPD to investigate the behavior of a binary mixture of A and B, exhibiting a demixing in a bulk phase, confined in slit-like pores with walls modified by the stripes of tethered brush of chains made of beads A. The main issue was to determine possible morphologies that can be formed inside the pore depending on the geometrical parameters characterizing the system (the size of the pore and the width of the stripes), cf. figure 6. In order to describe the observed morphologies, these authors calculated several characteristics, such as the density and local temperature profiles, the radii of gyration for the attached polymers, and the minimum polymer-polymer distances in the direction parallel and perpendicular to the pore walls. The summary of findings was presented as a sketch of the diagram of morphologies.
The study of reference [160] has been extended to the case of different (so-called in-and out-of-phase) arrangements of the stripes at two surfaces and varying composition of the confined mixture [161]. For in-phase narrow stripes and narrow pores, the case is reduced to a quasi two-dimensional micro-phase separation within the mixture confined between the polymer-rich lamellae on each wall. On the other hand, for narrow pores and moderate stripe widths, the in-pillar separation is quasi one-dimensional. In the case of wider pores and very wide stripes, the pore interior splits into the regions showing quasi twodimensional and bulk-like behavior. A number of solvent-mediated morphologies are observed when the stripes are weakly separated, both along the wall and across the pore. With an increase of the concentration of A beads in the mixture, they bridge up the striped brush in either or in both directions. Certain effects of morphology switching are observed. On the application side, the concentration driven morphologies can be fixed permanently via crosslinking, if a mixture contains a short crosslinker compatible with the polymer. Fixed morphology could be used as a nanotemplate or as a microfluid channel [161].
An interesting work on the application of DPD simulation technique to describe brushes on nanopatterned surfaces was published by Patra and Linse [138]. In the model considered by them, polymer chains were randomly attached onto an infinite long stripe along the Y axis with finite width ∆ along the OX axis. The polymers were treated as jointed chains composed of M spherical subunits (beads) connected by bonds. Each bead corresponded to a small part of a real polymer, approximately one Kuhn length [173], and was described by a soft repulsive potential. The simulations were carried out for a range of widths of the stripe, ∆, for several grafting densities ρ C and for several numbers of beads, M .
In figure 7 we show examples of the local densities of the beads as functions of X (along the surface) and z (in the direction perpendicular to the surface). The calculations were carried out for several values of ∆ (all distances are measured in the unit of the bead diameter). The "empty" space is occupied by a "bad" solvent. One should note solvent-mediated switching between lamellar and "meander"-like morphology in the bottom row of images. Reprinted with permission from reference [161]. Copyright CMP.
As ∆ is increasing, the following trends in the system were observed: (i) an increase in the brush height h, (ii) an expansion of the region with high density close to the grafting surfaces, and (iii) an appearance of a considerable overshot of polymers outside the grafting region. The polymer density fell off rapidly when the perimeter of the brush was reached.
Formally, the brush profile over a stripe can be evaluated by using AFM. However, the AFM tip, when approaching the surface, deforms the brush and, therefore, the height obtained from experiments differs from that for undeformed brush. Typical shapes of brush profiles from simulations and AFM measurements, as estimated by Patra and Linse [138], are displayed in figure 8. In contrast to the simulated profile, the experimental one does not display a clear plateau in the center of the stripe even though the brush is 10 ÷ 100 times as wide as it is high. Patra and Linse also present a strict derivation for the maximal decrease of the brush height outside the center. Since a polymer brush is a soft material, it bends under the force of the AFM tip. This behavior was also obtained by MD simulations of homogeneously grafted polymer brushes [174].
The simulations by Patra and Linse have been carried out in the domain M ∆ 5M . Their results suggest that this domain contains the transition from narrow polymer brushes to the regime of wide brushes. An improved understanding of this transition calls for a reassessment of the experimental conditions of the AFT technique, but even more exciting task would be to investigate the experimentally preferable orientation and conformational changes of chains near an edge by using, for example, labeled polymers. Many applications of nanopatterned polymer brushes would benefit from such a deeper understanding of the molecular arrangement in narrow polymer brushes.  The location of the stripe is indicated by the filled black area. Reprinted with permission from reference [138]. Copyright ACS.

Atomistic and semi-atomistic modelling
The strength of mesoscopic simulation methods, discussed in the previous subsection, is focused around their capability to cover nano-scale structure and the dynamics of polymer brush and the (mixture of) fluids contained within the pores. However, the soft interaction potentials, used in such studies, are typically of heuristic form and their relation to the real chemical substances is still an open question. Rigorous coarse-graining methods are still to be developed [175][176][177][178]. Therefore, the atom-based simulations have a twofold purpose: (i) to be used on smaller scale systems to serve as a feed for the estimates of coarse-grained interaction (so-called multiscale approach) and (ii) to be an independent tool for more "chemically detailed" simulation studies. An extensive literature is available on the applications of (semi-)atomistic models, reviewed earlier in references [143,144,[179][180][181]. Here, we outline only some of the works that are of interest or have been omitted in previous reviews.
One of the major tasks of computer simulations is to verify the analytic predictions from the SCF theory for grafted polymer layers. The MC simulation using the bond-fluctuation model was employed in [183] where the quantities describing the equilibrium structure of the brush are derived from the SCF theory and compared with the MC data without free parameter. In most cases, the results are in agreement with the SCF predictions. The causes of discrepancies are also discussed. The MD simulation of polymeric brushes in solvents of varying quality is presented in reference [91] where the resulting equilibrium structures are compared with the predictions of SCF approaches. The scaling of the brush height was found to agree with SCF predictions at temperatures T > T θ , where T θ is the θ temperature [2] of the solvent. For T < T θ , the agreement was poor, which was connected with the separation of the brush into solvent-rich and solvent-poor phases. This phase separation was found to disappear with an increasing grafting density and chain length. The effect of a weak attractive surface interaction on the structure of the brush was also studied.
The entropy of polymer brushes in a good solvent limit was evaluated numerically in reference [184]. The cases of a free polymer brush, a compressed polymer brush and two polymer brushes facing each other were studied. The total number of configurations, entropy, and the force acting on the polymer were calculated for various chain lengths and grafting densities by means of MC simulations using an efficient enrichment algorithm on a simple cubic lattice. It was demonstrated that the effect of an excluded volume plays an essential role and that the entropy of a free polymer brush drops as (spacing) −2.4 , while the entropy of a compressed polymer brush and two polymer brushes agrees very well with the analytical SCF theory by Milner, Witten, and Cates [185]. The analysis has been also extended to the case of a polymer solvent.
A MD simulation of polymeric brushes immersed in a melt of mobile polymer chains was presented [186]. The brush height and segment density profiles were evaluated for chains of length M grafted at one end to a solid surface, immersed in a melt of polymers of the same type. As the chain length P of free chains increases, there is a crossover from a wet to a dry brush in agreement with scaling and SCF theories. Since the interaction between all segments is identical, this crossover is driven purely by entropy. The simulation results were compared with earlier simulations. The exchange of chains between the brush and the bulk solution is of obvious relevance in understanding the thermodynamics and the kinetics of the formation of a grafted layer. In reference [187], by including a finite head-group adsorption energy, the bond-fluctuation model has been used to simulate

12601-23
the polymer brush in a good solvent with the effect of chain exchange between the brush and the bulk. The self-adjusted surface coverage was measured for different values of chain lengths and head-group energies. It was found that the polymer chains in a grafted layer are replaced by introducing shorter chains of identical head groups, which supports some experimental studies. The irreversible adsorption of single chains grafted with one end to the surface was also studied by combining the scaling theory with computer simulations in reference [188]. Using scaling arguments, the authors derived relations between the overall time of adsorption the characteristic time of adsorption, and the chain length. To support the analysis, the MC simulations were performed using the bond fluctuation model. The sequence of adsorption events is reproduced very well by simulations, and an analysis of various density profiles supports the theoretical model.
To clarify the differences in the properties of grafted chains at various grafting densities, a number of simulations were undertaken. In particular, the MC simulations were performed using bonded Lennard-Jones spheres. These simulations were aimed at observing the mushroom-to-brush transition at low grafting densities as well as the properties of a high density brush [189]. Through approximate scaling arguments, the location of the mushroom-to-brush transition was estimated and found to agree closely with the grafting density at which no more chains with mushroom dimensions can be inserted without overlap. Comparisons with analytical SCF theory for poor solvents were found to be much less favorable. In this case, the density profiles substantially differ from the SCF picture due to high segment densities and cluster formation. In reference [190], the static and dynamic properties of polymer brushes of moderate and high grafting densities were studied using MD simulations and the Langevin-type model. The fluctuation dynamics in the direction lateral to the surface was found to be well described by Rouse scaling. The dependence of the dynamics on the grafting density was well described by the dynamic behavior of thermal blobs for fluctuations in the direction both perpendicular and parallel to the surface. The pulling forces acting on the polymer bonds were also studied.
The knowledge of the dependence of the bond forces on the grafting density is important to develop strategies for making brushes with high grafting densities. The MD method allowed the authors of reference [190] to calculate the average vertical forces pulling the chains and the corresponding grafting energies. The binding potential between the nearest chains was the FENE potential, for details see reference [190]. The obtained pulling force as a function of the average distance from the surface for chains of M = 64 segments is displayed in figure 9. The results were calculated at different grafting densities (shown in the plot). At the top of the brush, the stretching forces vanish, as long as the grafting density remains moderate. At high grafting densities, even the end bond displays some residual stretching as a result of strong interactions with neighboring chains. Moreover, at high and moderate grafting densities, one observes a plateau which indicates the brush regime. Reprinted with permission from reference [190]. Copyright ACS.
The wetting behavior of polymers on top of a brush consisting of identical polymers has been studied using the computational methodology developed for a direct measurement of a wetting transition and its order via the effective interface potential. The method also permits to estimate the contact angles in the non-wet state and to study the adsorption isotherms [182]. In the absence of long-range forces, the system shows a sequence of non-wet, wet, and non-wet states as the brush density is increased. Having included attractive long-range interactions, we can make the polymer wet the brush at all grafting densities, and thus both first-and second-order wetting transitions are observed. The study highlights a rich wetting behavior when competing adsorbent-substrate interactions of different scales are tuned over a broad range.
Polyelectrolyte brushes have also been simulated using the MD method [191]. These were aimed at studying the fully charged polyelectrolyte brushes with salt, and the results are compared with SCF theory including finite stretching and volume effects [192]. The SCF approach was found to be capable of reproducing the brush heights at different grafting densities, salt concentrations and chain lengths at a semi-quantitative level. At high grafting densities, the density profiles obtained with both techniques exhibit a particularly close agreement, while at low densities, systematic deviations between their shapes are observed. The approximation of local electroneutrality, which the SCF approach is based on, was studied, and its implications were discussed.
It is instructive to compare the local net charges of the vertical layers of the brush (figure 10). This quantity serves for validation of the local electroneutrality assumption made in the SCF approach. Not surprisingly, this assumption is satisfied rather well when the Bjerrum length is large (i.e., for strong electrostatics) and less within the weak charge regime. Throughout the osmotic regime (we recall that in the osmotic regime, the concentration of counterions inside the polyelectrolyte layer is higher than in the bulk solution [193,194]), there exists a charge separation near the brush surface. The counterions form a charged layer just above the brush as a result of their osmotic pressure. It is this layer that generates the osmotic pull which contributes to the chain stretch. On the contrary, the local charge neutrality of the osmotic regime is rather well satisfied deep inside the brush. In the case of thick brush layers, when the surface contributes just a small fraction to the brush height, the approximation of local neutrality would turn even more accurate [191]. Reprinted with permission from reference [191]. Copyright ACS.
One of the topics extensively studied by using computer simulations is the adsorption of fluids on the surfaces modified by tethered chains [182,184,189] and the chromatographic process on columns filled with adsorbents modified by tethered oligomers [200][201][202][203][204][205][206][207]. A popular technique in the liquid chromatography is the Reverse Phase Liquid Chromatography (RPLC). The mobile phase in RPLC is a mixture of water and an organic solvent (usually methanol or acetonitrile). The stationary phase consists of a support material (usually beads of silica of micrometer size) and a bonded phase. The bonded phase is most often built of C18 alkyl chains. The retention of solutes occurs primarily within or on top of the bonded alkyl chains of the stationary phase.
The mechanism of retention in RPLC has been of interest for over three decades, but many of its aspects are still highly debated [208,209]. In general, solute retention in RPLC is driven by the free energy gradient on passing from the mobile phase into the stationary phase. One of the pioneering MD simulations of the transfer of a simple nonpolar solute from a water/methanol solvent mixture into a C18 stationary phase at room temperature was carried out by Klatte and Beck [195]. These authors, in addition to a detailed examination of the local environment of the solute, computed the excess chemical potential profiles. The free energy change was consistent in magnitude with that expected for hydrophobic transfer from water rich to oil phases, but specific interfacial effects drew into question the bulk partitioning models.
To elucidate the molecular-level structural features that control shape-selective separations, the MD was performed for chromatographic models with alkylsilane length and temperature analogous to actual materials [196]. Alkyl chain order was found to increase both with an increase of the chain length and with a reduction of the temperature, as observed experimentally. Chromatography models of various chain lengths and over a temperature range that represents highly shape-selective RPLC stationary phases contained a series of well-defined rigid cavities. The size and depth of these "slots" increased for the C30 models, which may promote the enhanced separations of larger size shape-constrained solutes, such as carotenoids [196].
According to modern theories of RPLC, two main questions regarding the retention mechanism appear: (i) is solute retention better described by partitioning into the stationary phase that resembles bulk liquid-liquid partitioning or by adsorption onto the stationary-phase surface? (ii) what is the driving force for the retention process, unfavorable interactions with the mobile phase or favorable interactions with the stationary phase? The answer of the above questions was considered in a series of papers by Siepmann and co-workers [201][202][203][204][205][206][207], who employed an advanced MC simulation method. In particular, in an interesting recent work [206], the effects of stationary phase and the solute chain length, as well as the solute adsorption and separation were studied. The surface was modified by grafting dimethyl triacontyl (C30), dimethyl octadecyl (C18), dimethyl octyl (C8), and trimethyl (C1) to silane surface. The modified surfaces and the bare silane surface were in contact with a water/methanol (33% of methanol) mobile phase that contained normal alkanes and alcohols. The simulations were carried out using the Gibbs ensemble and the details are given in [206,210,211].
The computer simulations were carried out at T = 323 K and pressure of 10 atm, i.e., at the conditions typical of RPLC experiments. Figure 11 (panel a) gives typical snapshots of the configurations of the molecules in the selected systems.
Simulating tethered molecules, Siempmann et al. [201][202][203][204][205][206][207] investigated how the structure of tethered layer changed with the length of the chains. In particular, they measured the dihedral angles, defined as the consecutive angles between carbons separated by two bonds and the normal to the silica substrate, as well as the so-called "gauche defects" [212] (defined as the fraction of dihedral angles in the entire chain deviating by more than 60 • from the angle of the trans-conformer. They also studied the local composition and stationary phase thickness. For this purpose, density profiles for water, methanol, and stationary phase carbon atoms were evaluated, cf. figure 11 (panel b).
The plots also show the Gibbs dividing surface (GDS, i.e., a plane defining the border between the mobile and stationary phases) and the width of the interfacial region. The latter was defined as the range of z-values where the total solvent density lies between 10% and 90% of its bulk value. However, due to very large oscillations in the total solvent density, it was not possible to determine a GDS and the interfacial region for the bare SiO 2 . One can see several similarities between the systems with the chains of different length. First, one observes a peak of the methanol local density just above the location of the GDS, i.e., there is a significant enrichment of the organic modifier in this region that extends deep into the interior of the tethered chains. The enrichment of methanol at the surface was similar for all lengths of chains. Second, there is a depletion of the total system density for all systems (see cyan lines), i.e., a dewetting effect at the surface.
However, there are also some differences between the system with different chains. For the C1 system the interface is very sharp and the methanol enrichment is very pronounced. The solvent density shows an oscillating behavior while for other chains it decays smoothly in the interfacial region. For C8  [206]. The stationary phase is shown as tubes with silicon in yellow, oxygen in orange, hydrogen in white and carbon in black. The mobile phase is shown in the ball and stick representation with oxygen in red, hydrogen in white, and carbon in blue. The analyte molecules are depicted by large spheres with oxygen in red, hydrogen in white, and carbon in cyan. (b): System composition as a function of distance from the silica surface. Stationary phase carbon density is shown in black, methanol in blue, and water in red. The total system and solvent densities are shown in cyan and purple, respectively. The GDS is shown by the orange dashed line and the interfacial region is shaded gray. Reprinted with permission from reference [206]. Copyright Elsevier. system, the interfacial width is larger than for the other systems, and one observes a more significant solvent penetration in the C8 system, as compared to the C18 and C30 stationary phases. This effect was interpreted by examining the snapshots in figure 11. The C8 phase is incapable of completely covering the silica surface, and thus the filaments of solvent molecules bridging the substrate and the mobile-phase region prevail. In the case of C18 and C30 phases, they form a continuous layer with respect to the solvent penetration, and solvent bridges are extremely rare for C18 phase. In other words, the C8 phase forms a more heterogeneous surface layer. It should be stressed that the water and methanol densities in the interior region of the C30 phase are about one order of magnitude smaller than the corresponding densities for the C18 phase.
The simulations by Siepmann et al. indicated, that the retention of butane increases with an increasing chain length, which is in agreement with experimental data. The nonpolar alkane solute shows a mixed retention mechanism with contributions from adsorption in the interfacial region and partitioning into the stationary phase. For the C1 phase, retention is dominated by adsorption. This phase is too thin to allow for partitioning. However, the adsorption of butane on a "hard" surface of C1 is different from adsorption on a more flexible, or "soft" surfaces involving longer chains. Siepmann et al. have discussed not only the mechanism of retention of butane but also other molecules, and their calculations have clearly demonstrated that advanced computer simulations can be very helpful in explaining the mechanism of chromatographic separation.
The absorption (confinement) of free linear chains into a polymer brush was studied with respect to chain size M and the Flory-Huggins type compatibility parameter χ with the brush by means of MC simulations at various grafting densities using a bead-spring model. Different concentrations of free chains were examined. Contrary to the case of χ = 0, where all species are almost completely ejected by the polymer brush irrespective of their length M , it was found that for χ < 0, the absorbed amount undergoes a sharp crossover from weak to strong absorption, discriminating between oligomers, 1 M 8, and longer chains. For a moderately dense brush, the longer species, populate predominantly the deep inner part of the brush, whereas in a dense brush they penetrate into the "fluffy" tail of the dense brush only.
Gyration radius R g and end-to-end distance R e of absorbed chains thereby scale with length M as free polymers in the bulk. Three distinct regimes of penetration kinetics of free chains were found regarding the length M , in which the time of absorption grows with M at a different rate [197].
To summarize this part, one could just cite some simulation works that target technology-related aspects of polymer brushes. For instance, the authors of reference [198] remark that "via computer simulations, we demonstrate how a densely grafted layer of polymers, a brush, could be turned into an efficient switch through chemical modification of some of its end-segments. In this way, a surface coating with reversibly switchable properties can be constructed. We analyze the fundamental physical principle behind its function, a recently discovered surface instability, and demonstrate that the combination of a high grafting density, an inflated end-group size and a high degree of monodispersity are conditions for an optimal functionality of the switch." In this way, atomistic simulations of polymer brush involving materials are expected to increase their predictive power and help to create a variety of smart, switchable, and multifunctional surfaces and thin films. Applications of these versatile and multifunctional brush coatings are envisioned in many areas including fluid control, microfluidics, and thin film sensors as well as some objects of interest in molecular biology [199].

Summary
To summarize, in this text we have reviewed the application of DF approaches and computer simulation methods for a description of fluid-tethered polymer brushes interfaces. Several illustrations of important findings were given. The topic is very ample and, consequently, we were unable to describe it in full detail.
Albeit the success in application of the theories outlined in this review, they will benefit of further improvements along different lines. In particular, a better description of the model of tethered polymer chains is necessary. Thermodynamic perturbation theory of the first order involved in the majority of present studies does not permit to capture structural peculiarities of several systems of interest. DF approaches encounter difficulties in an accurate description of attractive interparticle and long-range electrostatic interactions. A comparison of theoretical predictions with computer simulation data is, therefore, highly desirable to avoid incorrect interpretation of some theoretical results. As for technical issues necessary to be solved in the nearest future, it is of interest to include statistical mechanical analogue of mass action law into the available theories. This would permit to explore the effects of chemical association more in detail. We are convinced that theoretical developments in the area would permit to find novel phenomena faster and with less cost compared to other methods.
Our review also presents selected applications of computer simulations in the studies of systems involving tethered chains. First, we briefly described some results obtained from DPD coarse-grained simulations and then we outlined the application of microscopic MC and MD simulations to explain the mechanism of separation in RPLC. With an increase in performance of computer facilities, the role of simulations will increase, and the simulation methods in many cases will replace not only theoretical approaches, but also some experiments. This is particularly true for the cases of complex systems and for out of equilibrium phenomena, such as fluid flows. One of possible directions toward the development of simulation techniques to describe brush-fluid interfaces would be the application of multiscale modelling using field-theoretic methodologies [213]. According to those methods, some fragments of polymers are treated on a fully atomistic (or even ab initio) level and the information is used for parametrization of coarse-grained potentials [214]. Such a technique has already been used to describe the formation of micelles and vesicles in amphiphilic systems [215]. However, when performing coarse-graining, one should always keep in mind that the coarse-grained systems are physically different from the initial, molecular systems and that the coarse-graining methods cannot be used as a black box and require a thorough cross-checking.