Studying polymer solutions with particle-based models linked to classical density functionals : conon-solvency

We demonstrate the potential of hybrid particle-based models, where interactions are introduced through functionals of local order parameters, in describing multicomponent polymer solutions. The link to a free-energy-like functional is advantageous for controlling the thermodynamics of the model. We focus on co-non-solvency – the collapse of polymer chains in dilute mixtures with two miscible good solvents, having different affinities towards the polymer. We employ a simple model where polymers and solvents are represented, respectively, by worm-like chains and single particles. Non-bonded interactions are captured by a polynomial which is third order in local densities and can, therefore, describe liquid–vapour coexistence. The parameterisation of the functional benefits from an elementary mean-field approximation to the statistical mechanics of the model. The model provides a framework for Monte Carlo simulations using a particle-to-mesh algorithm. Studies with conventional generic bead-spring and all-atom models have demonstrated that co-non-solvency is caused by preferential binding of the better solvent (termed cosolvent) with polymer. Hence, segmental loops bridged by cosolvent molecules are formed, initiating polymer collapse. The mesoscopic hybrid model differs conceptually from the conventional microscopic descriptions. Yet, it reproduces the same co-non-solvency mechanism supporting its universality. Films of adsorbed ternary solutions, showing co-non-solvency in the dilute regime, are considered at high concentrations. In this case, chains do not collapse. The properties of loops and tails of the adsorbed polymer agree with early theoretical predictions obtained for concentrated binary solutions.


Introduction
Many fundamental studies and technological applications of soft matter involve polymers mixed with several solvents.To understand theoretically complex interrelationships 1 between structure, processing, and properties of multicomponent polymer solutions, molecular simulations are commonly required.These interrelationships often involve length and time scales which can be addressed only using mesoscopic descriptions: molecular models where large groups of actual chemical units are mapped on a single effective particle.Drastic coarse-graining reduces significantly the degrees of freedom and softens the interactions, making them comparable to the thermal energy.Therefore the computational time remains tractable.The difficulty in defining the mesoscopic interactions is one of the downsides of drastic coarse-graining.Standard bottom-up strategies 2,3 require reference data from computationally expensive simulations based on detailed models.Moreover, performing reference simulations of polymer solutions with new and complex formulations, [4][5][6] often requires adjustments of existing atomistic force fields.Finally, bottom-up strategies suffer from state-point and representability dependent parameterisation.
Top-down coarse-graining of polymer solutions is an alternative strategy where mesoscopic models are developed postulatively, aiming to reproduce by construction a set of known properties.Based on this pre-existing (limited) knowledge, topdown modelling can interpret and predict other features. 3hermodynamic properties of fluids play a key role in chemical engineering, [7][8][9] are extensively documented, and can provide guidelines for top-down modelling.Constructing and improving force fields using thermodynamic information is a common practice in atomistic and moderately coarse-grained (CG) models; the NERD 10 and MARTINI 11 force fields are typical examples.
Top-down particle-based models linked to continuum freeenergy-based descriptions allow one to introduce systematically [12][13][14] thermodynamic information.These approaches describe the molecular architecture of different compounds similarly to standard CG models.However, non-bonded interactions are templated by functionals [15][16][17][18][19][20] of local order parameters.Structurally, these functionals are similar to free energies in classical density functional theory (DFT) and phase-field models.Yet, they are conceptually different because they depend on instantaneous order-parameter fields (and not average).The functionals are transformed into mesoscopic force fields by expressing the instantaneous local order parameters through the coordinates of the CG particles.The hybrid method simplifies the construction of particle-based models with desired thermodynamics, because the thermodynamic properties are already determined by the templating functional.Simple mean field (MF) estimates, Self Consistent Field (SCF) theory, and classical DFT facilitate the choice of a functional suitable for reproducing the target thermodynamic behaviour.Importantly, introducing interactions through functionals of fields of local order parameters facilitates coupling with other continuum techniques.For example, hybrid schemes combining [21][22][23] dynamic SCF theory (based on quadratic density functionals) with the Lattice Boltzmann method have been used to explore the effect of hydrodynamic interactions on the kinetics of self-assembly in block-copolymer-based solutions. 22,23lthough the popularity of hybrid models in modelling Soft Matter is increasing, [12][13][14]24 their potential in modelling polymer solutions has not yet been fully explored. Curently, hybrid models of polymer solutions 22,23,[25][26][27][28][29][30][31][32][33][34] are build on functionals which are polynomials of local densities of different compounds.Most studies focus on phenomena in the liquid phase (such as self-assembly) and retain only second-order terms, which are equivalent to simple Flory-Huggins interactions.Significantly less studies incorporate polynomials of higher order, which are necessary to cover a broader range of phase behaviour.27,35 For example, using third-order density polynomials is the minimum requirement to have coexisting liquid and vapour phases. The few stuies with higher order functionals usually consider the solvent implicitly.Namely, the ''vapour'' phase is actually an integrated-out phase 12,25,26,28,29 of liquid solvent at equilibrium with an explicitly described polymer-rich phase.A single solvent component and an actual vapour phase have been introduced only in a very limited number of studies to describe processes in films of polymer solutions, such as drying 27 and solvent annealing.32 Here we expand previous studies and investigate whether simple mesoscopic hybrid models can capture complex phenomena sensitive to the formulation of multicomponent polymer solutions.We focus on the phenomenon of co-non-solvency, which occurs when a polymer collapses in mixtures of two competing miscible good solvents.This phenomenon has been extensively studied in various experimental setups.[36][37][38][39][40][41][42][43][44] From the simulation side a rather limited number of studies has been performed to understand this puzzling phenomenon, including chemical-specific 40,45,46 and generic simulations.[47][48][49][50] We note that in some of these studies 48 at least one of the solvents was described implicitly, in a mean-field manner.Co-non-solvency has also been addressed in a number of works based on analytical theories.36,[49][50][51][52][53][54] In previous works, two of us have demonstrated that the collapse of the chain is driven by preferential binding of the better solvent in the pair (named in the following cosolvent) with polymer monomers.49,50 The aggregation of cosolvent molecules onto the chains facilitates the formation of enthalpically favourable segmental loops. The pocess leads to the formation of a compact globule.This observation of preferential binding is also closely complemented by NMR and AFM experiments.41,44 Co-non-solvency depends 49 on the local structure of the liquid (i.e.energy density within the solvation shell).In hybrid models soft CG units overlap, there are many interacting neighbours, and the liquid structure is captured only on mesoscopic level.18 Furthermore, local liquid packing is affected by the specific method used to calculate 13 the interactions, e.g.particle-to-mesh (PM) schemes 15,18,28,55,56 or density distributions.[15][16][17]19,20 Due to these simplifications, it is not known beforehand whether hybrid models reproduce co-non-solvency and, if so, which features.
First we demonstrate that co-non-solvency can be indeed reproduced by hybrid models constructed using basic thermodynamic reasoning.We intentionally implement a simple hybrid model of ternary solutions which (i) describes explicitly the three components (polymer, solvent, and cosolvent), (ii) is based on third-order density polynomials, and (iii) is combined with a rudimentary PM scheme for calculating interactions.The modelling approach captures qualitatively most of the known features of co-non-solvency.Subsequently, we use this model to describe the same ternary solutions, coated on a solid substrate in the concentrated regime.In these films, we explore the properties of adsorbed layer and compare them with analytical predictions known from simple binary solutions.

Modelling concept
The modelling approach can be formalised through the canonical configurational partition function: The index a indicates the component type and n a denotes the number of molecules of each component.The coordinates of the s-th monomer in the i-th molecule are defined by r i (s) and N a is the total number of monomers in a molecule of component a.In this work, the three components are defined by a = p (polymer), s (solvent), and c (cosolvent).The integral Ð dR is taken over the possible realisations of the coordinates of all monomers.The ansatz is general enough to allow the description of branch-like polymers and solvents which are composed of several monomers, i.e. are chain-like.The architecture of molecules of This journal is © The Royal Society of Chemistry 2018 component a is described by the effective Hamiltonian of bonded interactions, H b(a) .The effective Hamiltonian of non-bonded interactions, bH nb , is expressed through a functional, f, of local densities, ra (r), of monomers belonging to the different components.In our study, we assume in eqn (1) that the molecules of each component are comprised of monomers of a single chemical species, e.g.polymers are homopolymers.Generalising to molecules with chemically different monomers is straightforward.The functional f depends on instantaneous values of densities (indicated by the ''hat'').Essentially, ra (r) are operators transforming a configuration of particle coordinates, R, into instantaneous fields of collective variables.To define the model, in addition to H b(a) and bH nb , the density operators must be specified.

Simple model for ternary solutions
In this section we describe a simple hybrid model for ternary mixtures.Though generic, the model is parameterised such that some thermodynamic properties are qualitatively similar to actual materials.

Particle-based description
Polymers are described using the discrete worm-like chain (WLC) model, also termed 57 as ''elastically-jointed chain''.Each chain has N p ''monomers'' placed at the ends and junctions of the WLC or, equivalently, N p À 1 segments with fixed length b.The bonded interactions are defined through the Hamiltonian: u i (s) is a unit vector oriented along the s-th segment of the i-th WLC, and e controls the stiffness of the WLC.In some applications, the segments of the discrete WLC can represent groups of microscopic degrees of freedom of an actual polymer.However, here b and N p bear no physical significance.Therefore, the discrete WLC encapsulates only two physicallyrelevant length scales in the ideal-chain limit: the Kuhn and the contour lengths.In the remaining part of the paper, we consider the Kuhn length l k (in the ideal-chain limit) as the fundamental length scale and express all lengths in terms of this quantity.The amount of WLC segments contained in one Kuhn segment, l k /b, is related to e via: This equation follows from a well-known relationship 57 between the persistence length l p and e, considering that l k = 2l p when e c 1.
Solvent and cosolvent molecules are represented by single beads, so that N s and N c equal unity.Therefore there are no bonded interactions associated with these two species.

Non-bonded interactions
We define non-bonded interactions through a third-order polynomial of local densities: The coefficients n ab and o abc are invariant to index permutation.The third-order density dependence can describe liquid/ vapor coexistence and is, therefore, essential for studying polymer films with a free surface.At a first glance, f in eqn (4) cannot penalise the formation of interfaces because it is a local functional.In fact, locality of f is not an issue for studying films of polymer melts or concentrated solutions: the explicit description of polymer connectivity via the WLC model, creates surface tension through the loss of conformational entropy at the liquid/vapor interface.This entropy loss is equivalent 58 to augmenting f by a Lifshitz term, 59 which depends on the square-gradient of polymer density.The simple simulation method in our study does not allow the implementation of local f when liquid/vapor interfaces can form for pure solvent or cosolvent (which are monomeric liquids).More clarifications follow in the next section.When required, interactions with the substrate are described using the generic potential: 60 where, z is the distance from the substrate.Here the interacting particle can be a WLC monomer, solvent, or cosolvent molecule.
The parameters L and l determine, respectively, the strength and the range of bU(z).

Simulation method
The hybrid model provides a framework for Monte Carlo (MC) simulations.We use the canonical ensemble and place the solutions into a simulation box with edges L g , g = x, y, z.The length of L g as well as the amount of polymer, solvent, and cosolvent molecules in the box, are specific to the studied system and will be discussed during the presentation of results in Section 4. Bulk solutions are modelled with periodic boundary conditions (PBC) in all three directions.When modelling films, there are no PBC along z-direction: the box is sealed at z = 0 and L z with the potential bU(z) and an enthalpically neutral hard wall, respectively.A MC algorithm requires the calculation of energies from particle coordinates.Calculating the bonded energy and the energy from interactions with the substrate is straightforward by direct substitution of coordinates into eqn (2) and (5).However, the calculation of non-bonded interactions requires instantaneous fields of densities.We obtain these densities from the coordinates of the particles, distributed in continuum space, through a PM scheme.The simulation box is discretised through a cubic lattice with cell-size DL.The densities ra (c) at a grid-node, with coordinates c, are calculated from: The function P(r i (s),c) assigns the monomers to the nodes of the grid.Here to obtain a simple model, we employ a PM scheme where each particle contributes only to a single lattice point.Namely, P = 1 if |r À c| g r DL/2, (g = x, y, z) and P = 0 otherwise.PM algorithms with higher order assignment functions are available, 13,28,55,61 where particles contribute to the density field at several grid points.On the one hand, such schemes are computationally more expensive.On the other hand, they effectively make the functional non-local 27 and allow, therefore, the description of liquid/vapour interfaces in dilute solutions.In addition, higher order PM approaches mitigate lattice-induced artefacts.
The integral defining H nb is discretised into a sum over the N cell nodes of the lattice, substituting Ð dr ¼ P N cell m¼1 DL 3 .After this step, MC sampling becomes straightforward.To sample the conformations of WLC we employ two standard algorithms: the slithering snake 62,63 (reptation) and the crankshaft 64 move, termed sometimes 65 ''flip''.During the reptation move, a chain is randomly selected.One segment is removed from a randomly chosen end of the chain and attached to the opposite end.We use the simplest version of the reptation algorithm, where the new orientation of the segment is selected randomly without bias.Since only one segment is involved, each reptation move can propagate the chain by maximum one bond length b.
The rotation angle 65 in the crankshaft (flip) move is chosen randomly from the interval [Àp,p].For solvent and cosolvent molecules we use small random displacements: each of the three Cartesian coordinates of the displacement vector is randomly chosen from the interval [ÀDl,Dl], where Dl = l k /5.Since these displacements were sufficient for the purposes of our study, we did not attempt to optimise sampling, e.g. by considering larger Dl.The new density r new a (c) in the cells affected by the MC move is calculated considering the coordinates of the displaced particles of species a at their proposed and old positions.Taking into account, the old density r old a (c) in the same cells, the difference in the nonbonded energy DH nb is obtained.Depending on the particle species changes of energy due to bonded interactions and substrate potential (when present) are calculated, as well.After all changes in energy are available, the Metropolis acceptance criterion is applied.In films (no PBC along z direction) the MC move is rejected straightaway if the displaced monomer, solvent, or cosolvent particles are found outside the boundaries of the box.
The lattice spacing sets the range of interactions and is a key parameter in PM schemes. 18,61We choose DL comparable to the length of the Kuhn segment, which is the fundamental length scale of our model, setting DL = 1.05.

Parameterisation strategy
To define bH b(p) we set the stiffness of the WLC to e = 2.99 which, according to eqn (3), leads to l k /b C 5. This empirical choice ensures that the number of particles per lattice cell in the liquid phase is on the order of 20 (cf.data provided further in the text).Because of this rather large number of interacting neighbours, the behaviour of the model in the liquid phase can be qualitatively estimated through elementary MF and the coefficients of the polynomial f follow from simple physical arguments.
Rudimentary MF approximation to the full statisticalmechanical description of the hybrid model given by eqn (1), leads to a Helmholtz free energy of homogeneous phases: The composition of the phase is described by the densities of polymer monomers, r p , solvent r s , and cosolvent r c particles.
In the following we use special notation to indicate the composition of liquid and vapour phases: r a(l) and r a(v) .The density of liquid phases of pure components will be indicated by r a(l) .The first term in eqn ( 7) is the entropic part, where Z a represents the single-molecule partition function of each component.Similarly to the standard Flory-Huggins model, this MF treatment assumes that Z a do not depend on composition.Hence, only translational entropy is present.From the simple Helmholtz free energy one can estimate observables such as pressure, P = ÀqF/qV, chemical potentials, m a = qF/qn a , and response functions, e.g.isothermal compressibility, 1/k T = ÀVqP/qV.
To parameterise f we use the observables: The thermodynamic observables for pure components or binary mixtures are obtained from these equations as a special case, setting the densities of the missing components to zero.We first determine n aa and o aaa controlling the thermodynamics of pure components.For this purpose, we consider liquid/vapor phase equilibrium in pure polymer, solvent, and cosolvent, fixing a priori the density and the compressibility of the liquid phase.We postulate that the density of WLC monomers r p(l) in the pure polymer liquid should be representative of classical polymer melts.It is convenient to express the monomer density as r p(l) C (l k /b)n o , where n o is the number of Kuhn segments found in the characteristic volume l k 3 .The quantity n o is a well known material parameter, which is related to rheological properties 66 and is often written equivalently as l k /p ( p is the packing length).For many polymers 66 1 r n o r 10 and we choose n o = 4.17, which corresponds to a standard polymer -atactic PMMA.For l k /b = 5, the number density of monomers in the polymer liquid is r p(l) = 20.85.Because our study is generic, we assume that solvent and cosolvent have different interactions with the polymer but are identical with respect to all other properties.Moreover, this simple choice guarantees that solvent and cosolvent are miscible in the entire range of molar compositions, which is the typical for solvent mixtures where co-non-solvency is observed. 36,50For simplicity we also assume that the number density of particles in the liquid phases of pure solvent and pure cosolvent is the same as for the polymer.That is, we set r s(l) = r c(l) = r (l) , where r (l) = 20.85.Hereafter we set T = 400 K, which is representative of temperatures where polymers are found in molten state.
Liquids are characterised by low compressibilities, i.e. typically k T B 10 À9 -10 À10 Pa À1 .However, hybrid models with small compressibility pose significant challenges to PM schemes.The intrinsic width of the interface of a polymer liquid (melt or solution) with its vapor is on the order of the Edwards correlation length x.In melts the Edwards correlation length [67][68][69] x m , is given by Here x p is the volume fraction of polymer segments.Using these expressions it is straightforward to estimate that k T B 10 À9 Pa À1 leads to x B 10 À2 Â l k .Therefore, realistic compressibilities lead to interfaces which are significantly narrower than DL and lattice-related artefacts are expected.To broaden the interfaces, we use for liquid polymer, solvent, and cosolvent k T = 10 À8 Pa À1 .Increasing the compressibility by an order of magnitude mitigates lattice-related artefacts, while the density contrast between vapor and liquid phases remains high (cf.discussion in Section 3.5).Using r p(l) and k T as an input, n aa and o aaa are estimated from the single-component limit of eqn ( 8)- (10).The calculations for the polymer are simplified considering that the pressure in polymer melts at equilibrium with their vapor is negligibly small and can be therefore approximated by zero. 70Substituting P = 0 and r p = r p(l) , eqn ( 8) and ( 10) are solved with respect to n pp and o ppp .For the solvent and the cosolvent the vapor pressure is not known and the complete condition of phase equilibrium is required, given by P( r . Unknowns in these equations are n aa , o aaa , and the saturation density of the vapor phase r sat a(v) (here a = s or c).The parameters n aa and o aaa obtained from the solution of the equations are listed in Table 1, while r sat a(v) C 0.5.Next we specify the parameters determining polymersolvent and polymer-cosolvent interactions: n pa , o paa , and o ppa .The coefficient n pa plays a key role in controlling the affinity of the solvent (cosolvent) to the polymer.The value n ps = À0.36508used to reproduce good solvent conditions is included in Table 1.This choice is based on trial simulations exploring the parameter space around n pa = (n pp + n aa )/2.At this reference point, the Flory-Huggins parameter w, to a first approximation, 12 equals zero, i.e. good solvent conditions are guaranteed.We note that this approximation neglects effects of composition on w from the third-order density terms in f.A more rigorous treatment of composition effects on w can be found elsewhere. 32,70In general, the isothermal compressibility k T obtained via eqn (10) will depend on composition.We reduce the dependency of k T on composition by suitable choice of o pss and o pps .Details are provided in the Appendix.
Because the solubility of polymer in cosolvent is higher than in solvent, n pc o n ps .The magnitude of n pc required to observe co-non-solvency can be estimated 45,49 from ternary mixtures of disconnected polymer monomers, solvent, and cosolvent particles.At infinite dilution, r p -0, the excess chemical potential of monomers follows from eqn (9) as a function of composition of the mixture: The composition is described via the volume fraction of the cosolvent x c = r c / r (l) .Based on eqn (11) we define a difference of excess chemical potentials Dm p(ex) (x c ) = m p(ex) (x c ) À m p(ex) (0).Observations in conventional particle-based simulations empirically suggest 45,49 that co-non-solvency occurs first when Dm p(ex) (1) is somewhere between À2 and À4kT.Setting n pc = À0.84 leads to Dm p(ex) (1) C À6.8kT which, as will be demonstrated in Section 4.1, reproduces co-non-solvency.Note that for each choice of n pc the parameters o pcc and o ppc are determined as in the case of solvent, i.e. they cancel the dependency of k T on composition.Table 1 lists these parameters for n pc = À0.84.
Because solvent and cosolvent particles are identical to each other (except from their interactions with the polymer) we set o ssc = o scc = o sss = o ccc .The single parameter related to simultaneous interactions of the three different species is determined via the combination rule o psc ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi o pss o pcc p .Regarding the interactions with the substrate, we set for solvent and cosolvent particles L = 0.For the interactions between the monomers and the substrate we employ L = 0.32 and l = 0.6.With this choice, the adsorption energy per unit area is estimated (roughly) to E ads =kT ' r pðlÞ Ð bUðzÞdz ' À5l k À2 .For a typical Kuhn segment length l k = 1.5 nm, this estimate leads to E ads /kT C À2.1 nm À2 which is comparable to magnitudes reported for actual materials. 60,71,72

Limitations of the generic hybrid model
We expect straightaway that the hybrid model introduced in the previous sections has a number of limitations.While the generic third-order density polynomial in eqn ( 4) is sufficient for our qualitative study, it is too simple to reproduce accurately the thermodynamics of most real systems, e.g.their equations-of-state.Such limitations have been discussed, for example, when modelling 35 alkane oligomers mixed with CO 2 or water on corrugated substrates. 73Various more complex free-energy models are available and are a basic tool for describing the thermodynamics of actual materials. 74pplying density functionals based on these models in hybrid modelling has significant potential but remains largely unexplored.For instance, the Sanchez-Lacombe equation-of-state has been only recently applied in hybrid modelling of homopolymer melts. 75In the past, third-order density polynomials have been applied 70 to systems forming hydrogen bonds.In this case they are analogous to Flory-Huggins models with negative Flory-Huggins parameters.A fundamental restriction of this approach is that it cannot capture the saturation of hydrogen bonds. 76his problem can be addressed using association models that introduce hydrogen bonding as an explicit reaction between particles.][78] In this case, the determination of the free-energy part which is related to the association of the particles presents a difficult combinatorial problem.This complication is avoided in hybrid simulations because they operate with ensembles of multiple, correlated chains.Therefore, it is possible 79 to describe association introducing reversible bonds between particles, e.g. through appropriate MC moves, similarly to standard particle-based simulations of equilibrium polymers. 80he parameters of the model are determined to satisfy MF equations linking them to thermodynamic observables used as an input.MF-based parameterisation is straightforward but does not allow one to completely control the thermodynamic properties of the hybrid model.The reason is that hybrid simulations incorporate fluctuations and correlation effects which MF neglects.Therefore, for a given set of parameters, the thermodynamic properties in hybrid simulations usually differ from their MF counterparts.Such deviations will be reported in Section 4. Therefore, the parameters of our model are reported in Table 1 with high accuracy only to facilitate the reproducibility of calculations discussed in Section 3.4.In practice, to obtain a hybrid model with qualitatively predictable thermodynamics, parameters (as those presented in Table 1) can be truncated -the MF-based parameterisation is anyway approximate.A simple rule of a thumb is to truncate the parameters such that after back-insertion into MF equations they reproduce the target quantities, e.g.densities, with an accuracy 0.1%.
As has been discussed in Section 3.4, making the liquid more compressible (increasing k T ) mitigates artefacts such as pinning of interfaces to the lattice.While this ansatz is sufficient for our generic qualitative study, we expect that reproducing actual vapour/liquid equilibria (VLE) of polymer solutions will require models with smaller k T .To identify the interdependence of compressibility and VLE, let us assume that the densities of the liquid and vapor phases of a pure solvent r s(l) and r sat s(v) , that are at equilibrium with each other at given temperature, are experimentally known.The coefficients n ss and o sss , necessary to reproduce this VLE point with a thirdorder density functional, are obtained by solving the set of two linear equations P( r s(l) ) = P(r sat s(v) ) and m( r s(l) ) = m(r sat s(v) ) (in this case n ss and o sss are the unknowns, cf.eqn ( 8) and ( 9)).Substituting these coefficients into eqn (10) delivers, within the model of a third-order density functional, the relationship between compressibility and density contrast of liquid and vapor phases: where x = r sat s(v) / r s(l) .Based on eqn (12), we illustrate how realistic VLE data can lead to low k T considering acetone as an example.For this standard organic solvent, two representative VLE points obtained 81

Bulk solutions
Before modelling co-non-solvency in ternary systems, we verify that our model reproduces characteristic behaviour of polymer solutions in single good solvent.We start from the dilute regime performing canonical MC simulations in cubic boxes with edge length L, containing n s solvent particles and a single polymer chain with N p = 661 (equivalent to 132 Kuhn segments).To model dilute bulk solutions, the total density in the box, r 0 , must be sufficiently high to exclude coexisting liquid and vapor phases.A first choice is to set r 0 equal to the liquid density, r s(l) , used for the MF parameterisation of bH nb .However, as has been discussed in Section 3.5 the thermodynamics in the particle-based simulations differs from the approximate MF solution.Therefore setting r 0 = r s(l) leads to the formation of small vapor ''bubbles'' motivating us to work with somewhat higher densities, r 0 C 22.56.Fig. 1 presents the log-log plot of the polymer structure factor, S(q), calculated from simulations in the dilute regime.The intermediate q range follows the scaling S(q) B q À5/3 (blue line), as expected for polymer chains in good solvent (dilute regime).At larger q, S(q) B q À1 (red line), manifesting the rod-like behavior of the chain on length scales comparable to the Kuhn segment.For these simulations we set L = 23, which is more than three times larger than the radius of gyration, R g(0) C 6.5, of the polymer in the dilute solution.
This journal is © The Royal Society of Chemistry 2018 We increase gradually the concentration replacing solvent particles by polymer chains containing an equivalent amount of monomers.This procedure preserves the total number density fixed to r 0 C 22.56.The main panel of Fig. 2 presents the average gyration radius R g as a function of polymer volume fraction x p (the average monomer number density normalised by r 0 ).As expected for polymers in good solvent, R g becomes smaller as x p increases and converges to the value for the melt (x p = 1).We compare our results with the scaling law 84 R g /R g(0) = (x p /x p *) Àf in the inset of Fig. 2. The overlap concentration is defined as x p * = 3N p /4p(R g(0) ) 3 r 0 , and f = (2n eff À 1)/(6n eff À 2).The scaling exponent n eff = 0.56, obtained by fitting the data for concentrations 1 o x p /x p * o 10, is not that far (B5% smaller) from the theoretically expected n eff = 0.588.The boundaries between dilute, semidilute, and concentrated regimes are not sharp but rather smooth crossovers. 85A plausible explanation for the observed deviation from n eff = 0.588 is that the system in Fig. 2 is found mostly in a crossover region.In fact n eff = 0.56 is within the range 0.54 r n eff r 0.63 reported in recent Brownian dynamics simulations 85 targeting crossover behaviour.
To verify that the hybrid model captures co-non-solvency, we perform hybrid simulations of ternary solutions in the dilute regime.The polymer has N p = 661 monomers and the initial configurations of ternary solutions with different composition are generated by replacing in the dilute polymer/solvent mixtures (cf.beginning of this section) some of the solvent particles with cosolvent.These configurations serve as an input for MC simulations.Fig. 3a presents the R g as a function of x c for n pc = À0.84 and two examples of weaker interactions: n pc = À0.7 and À0.56.For n pc = À0.84, the behavior of R g manifests the formation of a collapsed globular state at x c C 0.01.The inset in Fig. 3a illustrates the transition more clearly, focusing on the region 0 r x c r 0.1.The contraction of the polymer is substantially less pronounced for n pc = À0.7 and becomes insignificant for n pc = À0.56.Fig. 3b presents the probability distribution r(R g ) for the three considered values of n pc at x c C 0.01.For n pc = À0.84,r(R g ) is a narrow, strongly peaked function.In contrast, for n pc = À0.7 and À0.56, the distribution r(R g ) is broad manifesting a significant conformational freedom of the polymer chain.Note that the smallest R g for which r(R g ) 4 0 is located near R g C 2. This observation can be rationalised considering the properties of the polymer liquid Fig. 1 Single-chain structure factor, S(q), for a WLC polymer with N p = 661 monomers in a dilute mixture with solvent.The universal scaling laws S(q) B q À5/3 and q À1 expected for wavevectors probing (respectively) the length scales of the self-similar structure of the coil and the Kuhn segment are shown with blue and red lines.The mechanism of co-non-solvency in the hybrid model is the same as in the conventional all-atom 45 and generic 49,50 models.In all cases, the collapse of the chain is dictated by preferential binding of the cosolvent with the monomers of the polymer.Due to this binding, cosolvent particles form bridging contacts with monomers located far from each other along the contour of the same chain and initiate collapse.Fig. 4a presents a configuration of a collapsed chain (gray) together with cosolvent particles (red) in the regime of co-non-solvency (to improve visibility the solvent is not shown).To make evident that cosolvent segregates into the polymer-rich region, Fig. 4b reproduces the same snapshot showing only cosolvent particles.We emphasise that the collapse of the chain is not induced by the change of interactions between the components of the mixture as a function of composition (equivalent to a composition dependent Flory-Huggins parameter).As an illustration, Fig. 3c presents bDm p(ex) (x c ) for n pc = À0.56,À0.7, and À0.84 (black, red, and green lines).The plots demonstrate 50 that good-solvent conditions apply to the entire composition range; in fact the affinity of the solvent/cosolvent mixture for the monomers increases as cosolvent is added.
Fig. 5 demonstrates the normalised average radius of gyration 49,50 % R g = R g /R g (x c = 0) of the chain as a function of x c for n pc = À0.84 and several chain lengths, equivalent to N k = 10, 33, and 132 Kuhn segments.We observe that the range of concentrations where chain reopening occurs broadens as the length of the polymer increases.This observation follows the behaviour reported in conventional particle-based models of co-non-solvency 49,50 and has been explained through the larger flexibility of longer polymers to create loops. 49,50nterestingly, in standard particle-based simulations 45,49 of systems showing co-non-solvency, such as poly(N-isopropylacrylamide) (PNIPAM) in aqueous alcohol mixtures, the chain remains collapsed at a broader range of concentrations, x c , comparing to the hybrid model.This quantitative difference may originate from the longer range of interactions in the latter, defined by the lattice spacing.Generally, it is not easy to conclude how many cosolvent molecules participate in a bridging event in the PNIPAM polymer.However, we expect that in the hybrid model the contribution of a cosolvent molecule to bridging is more pronounced because it interacts with all polymer monomers located in the same lattice cell.Therefore the enthalpic gain from bridging contacts saturates in the hybrid model at lower x c .The polymer coil then swells again to increase conformational entropy.We observe in Fig. 3a that in hybrid simulations the window of concentrations, where the chain is collapsed, becomes more narrow as the contrast between solvent and cosolvent is reduced (n pc is closer to n ps = À0.36508).This observation matches the trends reported in standard particle-based simulations 50 and experiments. 36,37o verify that our simulations of co-non-solvency are free from pronounced mesh-related artefacts, we have calculated the average density of the polymer chain with respect to its center-of-mass (COM) in the expanded, x c = 0, and collapsed, x c = 0.01, states.A two-dimensional contour plot of the average density on a plane intersecting the COM is presented for these two states in Fig. 6a and b, respectively.The contour plot for the expanded chain is noisy due to strong fluctuations.In contrast, as has been already demonstrated in Fig. 3b, the fluctuations of chain conformations are significantly suppressed in the collapsed state, leading to a smooth distribution.In both cases  the isodensity lines are approximately circular, indicating that there is no strong influence of the grid on the interfaces.
The structure of the collapsed chain is examined in more detail in Fig. 7, presenting the average profiles of the densities of the three components, r p (r) (polymer), r s (r) (solvent), and r c (r) (cosolvent), as well as the total density of the system, r t (r).The profiles are calculated as functions of distance r from the COM of the WLC.Due to the strong preference of cosolvent for the polymer and increased compressibility of the model, we observe that in the vicinity of COM the total density rises by B20%.In agreement with Fig. 4, there is a significant increase in the density of cosolvent particles inside the collapsed coil.In simulations of co-non-solvency with standard particle-based models, the interface between the collapsed polymer and the surrounding solution is sharp. 45,49However, r p (r) in Fig. 7 has a rather broad decay, comparable to 2DL.We believe that this broadening can be explained considering the lattice-based definition of interactions in the hybrid model.Within this approach, monomers are free to move within a single lattice cell.This translation contributes a scale BDL to interface broadening.In fact, some registration of monomers within a lattice cell is expected 13,61 due to self-interactions, but is not sufficiently strong to cancel translational freedom.
Observing the same mechanism of co-non-solvency in the hybrid model and in standard particle-based models supports strongly its generic nature.Here we employ a hybrid model based on a third-order functional but we expect that using quadratic density functionals should be sufficient to observe co-non-solvency.Hybrid models with quadratic functionals are equivalent to using short-ranged pairwise potentials. 18In this scope, our expectation is corroborated by early simulations 47 where co-non-solvency was captured using a lattice model of polymer solutions with simple nearest-neighbor interactions.The advantage of having a third-order density polynomial is that the model can be directly transferred to studies of polymer films exposed to solvent vapour.

Films of binary and ternary solutions
In this section, we transfer the hybrid model parameterised to describe co-non-solvency in bulk dilute solutions (n ps = À0.36508 and n pc = À0.84) to study films in the concentrated regime.Unlike bulk solutions, the amount of polymer, solvent, and cosolvent molecules in the simulation box is now chosen to allow coexistence of liquid and vapor phases.Rectangular boxes are employed where the lateral dimensions depend on the length of the considered polymer chains; for N p = 166 we use L x = L y C 11 while for N p = 661 we set L x = L y C 23.These choices ensure that the lateral dimensions are approximately two times larger than R g(0) .In all cases, L z C 105.The initial configurations have a ''two-slab'' structure.The first slab is comprised of 400 stretched WLC (for all considered N p ) and is placed on the hard substrate.The second slab contains randomly mixed solvent and cosolvent particles and is placed on top of the polymer layer.The remaining part of the box is empty.The thickness of the two slabs is consistent with the density r (l) .During equilibration, these initial configurations evolve into adsorbed films of polymer solutions at equilibrium with solvent/cosolvent vapor.A typical equilibrated configuration is presented in Fig. 8a.For test cases initial configurations with different geometry were realised, to verify that the results are not affected by the initial setup.
It is instructive to compare the vapour/liquid equilibria (VLE) in simple MF treatment of our model and actual simulations.For polymers with N p = 166, we consider adsorbed films of binary polymer/solvent and polymer/cosolvent mixtures with different compositions.Fig. 8b shows a representative density profile orthogonal to the substrate for a polymer/solvent film.During simulations we have not observed evaporation of polymer molecules from the film, in agreement with Section 3.4 where the vapour pressure of the polymer was approximated by zero.Setting r p(v) C 0, the MF VLE for polymer/solvent mixtures are obtained solving the set of equations m s (r p(l) ,r s(l) ) = m s (r s(v) ) and P(r p(l) ,r s(l) ) = P(r s(v) ).The MF predictions for the total density, r t = r p(l) + r s(l) , and the polymer density, r p(l) , in the liquid phase as a function of r s(v) are presented in Fig. 9a with black solid and dashed lines, respectively.The MF VLE for polymer/cosolvent mixtures are calculated in a similar way and are presented in Fig. 9a with red solid and dashed lines.The densities r t and r p(l) in the bulk of the films as a function of r a(v) in simulations are shown with black (polymer/solvent) and red (polymer/cosolvent) symbols.The maximum dilution of the films that can be realised in our simulations is limited by the saturation density of the vapour phase, r sat a(v) .Above r sat a(v) solvent (cosolvent) droplets form, which cannot be described by the local functional of the current model.In simulations, the volume fraction of the polymer in the film at r s(v) C r sat s(v) is about 65% (for this estimate we use the approximate MF value r sat a(v) B 0.5).Films with cosolvent can be diluted to lower concentrations before r sat c(v) is reached, due to higher affinity of cosolvent towards polymer.However, for symmetry, we introduce the 65% limit in these simulations as well.
For similar composition of the liquid phase, i.e. comparable polymer concentration, the deviation between VLE in MF and simulations is smaller in polymer/cosolvent than in polymer/ solvent films.In other words, the deviation becomes smaller the better the solvent is.This observation suggests that the deviations in the VLE are largely due to the insufficiency of the MF approximation in describing the thermodynamics of the gas phase.To support this statement, we present in Fig. 9b the chemical potential bm(r s(v) ) of the solvent in the vapour phase as a function of r s(v) calculated from MF eqn (9) (line) and grand canonical MC simulations with the hybrid model (symbols).The deviations between r s(v) corresponding to the same chemical potential in MF and simulations increase as the gas phase becomes more dense.For small r s(v) the two descriptions are close to each other because the vapour phase is well approximated by the ideal-gas limit.Fig. 9b explains why the VLE in MF and simulations match each other better in polymer/cosolvent than in polymer/solvent films.To achieve low polymer concentrations in polymer/cosolvent films, small vapour pressures are sufficient.
In this case the vapour in both MF and simulations is close to an ideal gas.In contrast, the dilution of polymer/solvent films requires higher vapour pressures, where the ideal-gas limit is not applicable and MF fails.
Overall, Fig. 9a demonstrates that in polymer/solvent films, r t depends only weakly on composition.However, in polymer/ cosolvent films (where cosolvent interacts strongly with polymer) r t increases rapidly by about 14% at small concentrations of cosolvent in the film (corresponding to low r c(v) ).r t decreases slowly from this maximum value, as the film is diluted further.Exploratory MF calculations demonstrate that the change of r t with composition can be reduced for a model with realistic compressibility (data not shown).Of course, performing simulations with such a model requires more advanced techniques, comparing to the zero-order PM scheme in our study (cf.Section 3.5).
We consider ternary mixtures with different amounts of solvent and cosolvent to obtain films with various compositions and degrees of dilution.The example shown in the inset of Fig. 8b illustrates that the density profiles of these films are similar to the case of binary systems.Fig. 10 quantifies the dimensions of the polymer chains as a function of distance z of their COM from the substrate, in films of binary and ternary mixtures containing B65% polymer with N p = 166.Specifically, Fig. 10 presents the components of the root mean-square radius of gyration parallel, R 8 g (z), and perpendicular, R > g (z), to the substrate.Both quantities are normalised to the radius of gyration of the bulk melt.Although all films contain the same percentage of polymer, the films containing cosolvent are thinner.In these films, following Fig. 10, the polymer/vapor interface is located at z C 30 as opposed to z C 40 in the polymer/solvent film.This difference stems from the increase of the density in mixtures with cosolvent, as has been already discussed in connection with Fig. 9a.Polymers located close to the substrate and the polymer/vapor interface show an expected behaviour: 64,86 the normalised R > g (z) is significantly smaller than unity indicating ''flattening'' of chain conformations.In all three mixtures, the components of the radius of gyration indicate the same weak swelling of polymer chains in the bulk of the film, deviating from unity by about 5%.This behaviour is in agreement with the high polymer concentration in the films.Moreover, for films of ternary mixtures containing B65% polymer  (or higher) we observe that the chains remain swollen irrespective of the fraction of the cosolvent.Because polymer coils overlap strongly, the majority of enthalpically favourable bridges, created by cosolvent particles, will be between monomers belonging to different chains.Therefore, the entropically costly looping of chains is avoided and they do not collapse in the concentrated regime.
We quantify the structure of the film near the substrate using standard descriptors 87 related to the adsorbed polymer layer.The later is formed by segments belonging to adsorbed chains; a WLC is defined as adsorbed, if the distance of at least one monomer from the surface is smaller than DL.Density profiles of segments belonging to loops, r loop (z), and tails, r tails (z), of adsorbed WLC are presented in Fig. 11a and b.The plots consider a few representative films, for N p = 166 and 661, aiming only to illustrate the significant effect of composition and chain length on the profiles.The main result appears in Fig. 12a and b, presenting r loop (% z)/r p and r train (% z)/r p calculated in binary and ternary systems, for both N p , in the entire range of concentrations and compositions addressed in our simulations.Here r p is the density of polymer monomers in the bulk of the film and % z = z/R g(m) , where R g(m) is the average radius of gyration in the pure melt.This normalisation is motivated by the analytical theory of Semenov et al. 88 predicting that in the concentrated regime and z c R g(m) , the density profiles of loops and tails are described by the universal dependencies: Indeed the collective plots in Fig. 12a and b reproduce a mastercurve behaviour.The tail of the master curves can be described by eqn (13) (dashed lines) without any fitting parameter -we substitute for R g(m) the gyration radius which was obtained from simulations of melts.Fig. 13a-c present the probability distributions for the length of the loops, P loop (n k ), tails, P tails (n k ), and trains, P trains (n k ), calculated in binary and ternary systems with N p = 166.The lengths are expressed in terms of the equivalent number of Kuhn segments n k .We observe that in the concentrated regime the dilution has almost no effect on the statistics of conformations of adsorbed polymers.Motivated by Fig. 10, we compare the behaviour of probability distributions for loops and tails for large n k with laws known from adsorbed homopolymer melts: 89,90 The pre-factors C 1,2 are considered as fitting parameters and we note that the term BN k À1 n k À1/2 in P loop (n k ) takes into account the finite length of the chains.The theoretical predictions are shown in Fig. 13a and b by dashed lines and are found in good agreement with the simulation data.

Conclusions and outlook
We have demonstrated that a hybrid model, based on a simple density functional and a rudimentary (zero-order 13 ) PM scheme can reproduce the intriguing phenomenon of co-non-solvency.The local structure of the liquid in the hybrid model and in standard microscopic models -where co-non-solvency has been observed so far -differ significantly.The differences stem from drastic coarsegraining, the PM calculation of non-bonded interactions, and additional simplifications, such as higher compressibility.Nevertheless, the hybrid model reproduces the basic mechanism of co-nonsolvency proposed 45,49,50 in standard simulations.We observe that the collapse of a polymer chain in a solvent/cosolvent mixture, at small concentrations of cosolvent, is induced by preferential binding of cosolvent with polymer monomers.The binding promotes chain loops to maximise favourable contacts between polymer and cosolvent.Observing the same mechanism of co-non-solvency in the hybrid model and in standard particle-based descriptions supports strongly its universality.The contrast in the affinity of solvent and cosolvent towards the polymer required for co-non-solvency in the hybrid model, is comparable to conventional simulations.In hybrid simulations the chain remains in a collapsed state at a smaller range of compositions (than in conventional simulations).We also observe that, for the considered chain lengths, the interface between the collapsed polymer and the surrounding solvent is broader.We argue that the lattice-based definition of interactions contributes substantially to both effects.The ternary solutions, showing co-non-solvency in the dilute regime, were also considered at high concentrations.Benefiting from the fact that the polynomial functional underlying the hybrid model captures liquid-vapor coexistence, 27,35 the concentrated solutions were studied as films adsorbed on solid substrates.This geometry is relevant for coating applications.For polymer concentrations which are at least an order of magnitude higher than the overlap concentration, x*, we do not observe co-non-solvency.Due to substantial overlap of polymer chains, cosolvent molecules create enthalpically favourable bridges mainly between intermolecular monomers.Therefore, the collapse of individual polymer molecules is avoided.We quantified the structure of the adsorbed layer through standard descriptors. 87The properties of loops and tails in our ternary solutions match theoretical predictions obtained for concentrated binary solutions 88 and melts. 89,90verall, our study illustrates the potential of simple hybrid models in studying complex phenomena in multicomponent polymer solutions on qualitative level.However, the development of hybrid models reproducing quantitatively thermodynamic data of actual polymer solutions requires more complex approaches.For example, simple MF estimates of the thermodynamic behaviour of the hybrid model provide approximate guidelines for their parameterisation.In general, however, the thermodynamics of the hybrid model in particle-based simulations will be quantitatively different from the MF approximation.To reproduce quantitatively target thermodynamic properties, functionals that are based on more complex free-energy models 74,91 (rather than third-order density polynomials) and methods for systematic refinement of parameters estimated from MF are required.
from atomistic simulations (based on the OPLS force-field) are: (i) MW s r s(l) C 778 kg m À3 and MW s r sat s(v) C 1.46 kg m À3 at T = 300 K, and (ii) MW s r s(l) C 657 kg m À3 and MW s r sat s(v) C 18.2 kg m À3 at T = 400 K. MW s C 0.058 kg mol À1 stands for the molecular weight of the acetone (in contrast to the other parts of our paper, in the example for the acetone the number densities are expressed in mol m À3 for consistency with ref. 81).Eqn (12) delivers for these two VLE points k T C 4 Â 10 À9 and 9 Â 10 À9 Pa À1 , respectively.Liquids with k T B 10 À9 can be modelled 82 by combining density functionals with density distributions defined in continuum space, which avoids a PM discretisation.Multibody Dissipative Particle Dynamics 16,17,83 is an example of such strategies.

Fig. 2
Fig. 2 Main panel: average radius of gyration R g of WLC polymers with N p = 661 monomers in mixtures with solvent as a function of volume fraction of polymer monomers, x p .Inset: Reduced radius of gyration R g /R 0 g as function of reduced concentration x p /x p *.Here R 0 g is the radius of gyration in the dilute regime and x p * is the overlap concentration.The fit with the scaling law R g /R g(0) = (x p /x p *) Àf (see text for details) is shown with red dashed line.

Fig. 3
Fig. 3 (a) Average radius of gyration, R g , of a single WLC polymer (N p = 661) in a ternary mixture as function of cosolvent concentration x c .Three different strengths, n pc , of polymer/cosolvent interactions are considered; for all cases the polymer/solvent interaction is fixed to n ps = À0.36508.The region of x c where co-non-solvency occurs, is shown enlarged in the inset.(b) Probability distribution r(R g ) of the instantaneous radius of gyration at x c C 0.01 for the three cases of polymer/cosolvent interactions presented in (a).(c) Excess chemical potential of polymer monomers, bDm p(ex) , at infinite dilution, obtained as a function of x c from the MF eqn (9).The excess chemical potential at x c = 0 serves as a reference point.

Fig. 4
Fig. 4 (a) Configuration of a collapsed WLC (gray) with surrounding cosolvent particles (red) in the regime of co-non-solvency.To improve visibility, solvent particles are not shown.(b) The snapshot from (a) is reproduced showing only cosolvent particles, to make evident the segregation of cosolvent into the polymer-rich region.Fig. 5 Normalised radius of gyration % R g = R g /R g (x c = 0) as a function of cosolvent concentration x c for three different chain lengths expressed in number of Kuhn segments N k .For clarity, a small region of x c is enlarged in the inset.

Fig. 6
Fig. 6 Contour plots of the average polymer density on a plane intersecting the COM of a single WLC in dilute solutions.(a) Expanded state in binary solutions (equivalent to x c = 0).(b) Collapsed state in ternary solutions at x c = 0.01 (regime of co-non-solvency).

Fig. 7
Fig.7Density distributions, r a (r), as a function of distance from the COM of a collapsed chain in the state of co-non-solvency.Here a = t denotes the total density of the system, while the densities of polymer, solvent, and cosolvent are specified by a = p, s, and c.

Fig. 8
Fig. 8 (a) Configuration of a film of a binary solution, where polymer chains and solvent particles are shown with gray and blue colours, respectively.(b) Main panel: density profiles for a polymer/solvent (P + S) film containing 65% of polymer.Inset: Density profiles for a polymer/ solvent/cosolvent (P + S + C) film containing 65% of polymer.

Fig. 9
Fig. 9 (a) Comparison of vapour/liquid equilibria in polymer/solvent (black colour) and polymer/cosolvent (red colour) binary mixtures.MF predictions for the total density, r t , (solid lines) and the polymer density, r p , (dashed lines) in the liquid phase are presented as a function of vapour density of solvent, r s(v) , or cosolvent, r c(v) .The same quantities obtained from simulations are shown with squares and circles.(b) Chemical potential bm(r s(v) ) of solvent particles in the vapour phase obtained as a function of r s(v) from MF eqn (9) (line) and grand-canonical MC simulations (symbols).

Fig. 11 (
Fig. 11 (a) Density profiles, r loop (z), of segments belonging to loops of adsorbed chains for a few representative systems as indicated by the legends.The percentage indicated in the legends, refers to the polymer volume fraction in the film.(b) Same as (a) but for density profiles, r tail (z), of segments belonging to tails.

Fig. 12 (
Fig. 12 (a) Universal plot of renormalised density profiles, r loop (z)/r p , of segments belonging to loops of adsorbed chains as a function of rescaled distance, % z = z/R g(m) from the substrate.r p is the density of polymer monomers in the bulk of the film and R g(m) is the average radius of gyration in the pure WLC melt.The plot contains data representative of the chain lengths and compositions considered in our simulations.The percentage in the legends indicates the polymer volume fraction in the film.The analytical theoretical prediction from eqn (13) is shown with dashed line.(b) Same as (a) but for the case of tails.

Fig. 13
Fig. 13 Probability distributions for the length of (a) loops, P loop (n k ), (b) tails, P tails (n k ), and (c) trains, P trains (n k ), calculated in binary and ternary solutions containing WLC with N p = 166.The lengths are expressed in terms of equivalent number of Kuhn segments n k .The analytical theoretical predictions for large n k from eqn (14), are shown with dashed lines.

Table 1
Representative set of parameters for non-bonded interactions