Using field theory to construct hybrid particle-continuum simulation schemes with adaptive resolution for soft matter systems

We develop a multiscale hybrid scheme for simulations of soft condensed matter systems, which allows one to treat the system at the particle level in selected regions of space, and at the continuum level elsewhere. It is derived systematically from an underlying particle-based model by field theoretic methods. Particles in different representation regions can switch representations on the fly, controlled by a spatially varying tuning function. As a test case, the hybrid scheme is applied to simulate colloid-polymer composites with high resolution regions close to the colloids. The hybrid simulations are significantly faster than reference simulations of a pure particle-based model, and the results are in good agreement.


I. INTRODUCTION
Multiscale modeling is one of the central challenges in many areas of materials science [1][2][3]. The properties of modern materials are often determined by an interplay of structural features and processes on length scales that span several orders of magnitude. For example, many materials are heterogeneous on a nano-or micrometer scale and filled with "defects" -internal interfaces, droplets of a different phase, or nanoparticle fillers. Theoretical descriptions must account for the microscopic structure close to these defects as well as the larger scale structure of the "bulk" medium surrounding the defects [4]. To study such systems, multiscale modeling approaches have been developed and pursued for several decades, which employ a hierarchy of models to describe the material properties at different coarse graining levels [5]. One crucial issue in this context is the coupling between models. The traditional approach has been to couple them "vertically", i.e., simulations of different models are run independently and linked by parameter heritage. Nowadays, "horizontal coupling" schemes are attracting growing interest, where regions of different resolution coexist within one single simulation system [6]. In particular, the adaptive resolution models [7][8][9][10][11], which allow free diffusion of particles between regions of different resolution, are able to dynamically couple information and to account for density fluctuations and flow. The adaptive scheme is suitable for systems with small regions requiring detailed investigation, while the remaining large part only needs a computationally cheaper coarse-grained description. Such systems are ubiquitous in soft materials, e.g., chemical reaction systems, polymer solutions and melts with interfaces, or composite materials.
On the microscopic side, materials are typically represented by particle-based models (atomistic or coarsegrained). On the macroscopic side, continuum models are commonly used (elastic models, phase field models, hydrodynamic models). While horizontal coupling schemes have been developed both within the "particle world" and the "continuum world", linking the two still remains a challenge. Hybrid particle-continuum schemes have been proposed where certain molecules or components are treated permanently at the particle level, and others permanently at the field level [12][13][14]. Other examples of coupled schemes are "Single Chain in Mean Field" simulation methods, where particles move in the dynamically updated mean field of the surrounding particles [15][16][17], or 'heterogeneous multiscale' schemes where particle simulations are used to adjust the parameters of a continuum simulation on the fly [2,18]. However, apart from proposals for simple liquids [19][20][21], the present authors are not aware of a general scheme for complex fluids that would allow one to treat different regions of space at either particle or continuum level in an adaptive resolution sense.
With the present paper, we aim at closing this gap. We propose a method to generate adaptive resolution schemes that link particle and continuum representations of the same complex fluid in a formally exact manner. Together with existing adaptive particle-particle and continuum-continuum coupling schemes, our method could potentially pave the way to integrated multiscale treatments of complex fluids from the atomistic to the macroscopic scale.

II. BASIC CONCEPT OF THE APPROACH
Our starting point on the particle side are models of the Edwards type, which can be defined in terms of local densities. This implies, in particular, that the interaction potentials are soft, i.e., molecules can penetrate each other. Although the Edwards models were originally introduced in the context of analytical theory [22], they also proved to be efficient models for computer simulations [16,[23][24][25]. The partition function of an Edwardstype model can be rewritten exactly as a fluctuating field theory [26], either by applying a Hubbard-Stratonovich transformation (if the interactions are purely quadratic in the densities), or, more generally, by inserting unities (delta functions) in a Faddeev-Popov way [27][28][29][30]. Such fluctuating field models have also been studied by computer simulations with considerable success [31][32][33]. Even more importantly, fluctuating field theories lend themselves to mean-field approximations, thus providing a natural link between Edwards models and popular density functional theories for complex fluids such as the Self-Consistent Field (SCF) theory [27,30] or dynamic density functional theories [29,34,35]. These so-called "molecular field" theories are nonlocal continuum models, which can be used directly for mesoscale simulations of complex fluids [36,37], and which also provide an excellent starting point for systematic derivations of simpler phase field theories [38].
Thus every particle-based Edwards model has a continuum model partner, i.e., the corresponding molecular field model, which is equivalent apart from a meanfield approximation. We note that the mean-field approximation becomes accurate in dense systems, which is where the transition from a particle-based to a fieldbased model is most attractive. Moreover, the effect of fluctuations can often be included to some extent even in a molecular field simulation [29,32].
Our adaptive resolution scheme exploits this correspondence between Edwards models and molecular field models . We will construct a hybrid model that combines particle and field representations of the same molecules, and a simulation scheme to switch between representations depending on the position in space. The switching probability is controlled by a spatially varying virtual field ∆µ(r), which can be chosen at will. As an example, we will study a polymer-colloid composite, with ∆µ(r) chosen such that the particle representation dominates close to the colloids, and the field representation far from the colloids.
We will now describe the basic idea of our approach. Technical details are given in the Appendix A. For simplicity, we consider a one-component system of n polymers (labelled α) with N monomers (labelled j). Our starting point is the canonical partition function where the integrals run over all monomer positions R αj , H 0 denotes the Gaussian spring energy of the chains, and H nb describes the non-bonded interactions in terms of an Edwards Hamiltonian. Here and throughout, the energy unit is chosen 1/k B T ≡ 1. Based on the partition function (1), the hybrid particle-field model is now constructed in three steps.
In the first step, the polymer chains are partitioned into two different (virtual) species, which we name pchains and f-chains. This is done by attaching an additional virtual variable τ αj ∈ {0, 1} to each monomer. A chain α is said to be an f-chain if j τ αj = 0, otherwise it is called a p-chain. The virtual variables τ αj are introduced by inserting the exact identity 1 ταj =0 exp τ αj ∆µ(r) − ln e ∆µ(r) + 1 = 1 (2) in the partition function, Eq. (1). This couples them to the virtual field ∆µ(r), and the latter can be used to control the fraction of f-and p-chains at a given position r. We note that we are free to choose the field ∆µ(r) as we like, since the identity, Eqn. (2), is exact. The second step is to treat the p-chains and the fchains by different representations. We keep the particle description for the p-chains, but convert the description of f-chains into a field representation. This is done in the usual Faddeev Popov way by inserting appropriate identity operators (see Refs. [27][28][29] or Appendix A). As a result, the particle degrees of freedom of the f-chains are replaced by fluctuating fields φ f and ω f .
The resulting expression for the partition function is formally equivalent to Eq. (1), but it cannot be sampled efficiently. Therefore, the third step is to introduce convenient approximations that speed up the numerical calculations. Here we use a saddle point evaluation [27,35,39] of the ω f integral, and just keep the φ f fields. Such a mean-field type treatment only influences the contributions from the f-chains. It amounts to a kind of 'coarsegraining' in the low-resolution region. If the density of the medium there is high, the mean-field approximation is known to describe the system very well. The physics in which one is interested, however, is extracted from the high-resolution part where the polymers are still represented by particles, for which no approximations were used.

III. APPLICATION EXAMPLE: POLYMER-COLLOID-NANOCOMPOSITE
As an application of our hybrid model, we study a complex composite system containing two nanocolloids that are coated uniformly with homo-brush polymers and immersed in a melt of n t A-B diblock copolymers. Each free polymer consists of N = 20 monomer beads, with N A = 10 A-beads and N B = 10 B-beads, and each brush polymer contains 10 monomer beads. One colloid is coated with A-homopolymers, the other with Bhomopolymers. The non-bonded Edwards Hamiltonian for this system is given by where the Flory-Huggins parameter χN = 9 measures the incompatibility of monomers A and B, κN = 10 is the compressibility. The configuration dependent densities of monomers A and B are denoted byφ A andφ B , and φ 0 is the reference monomer density in the bulk fluid. Furthermore, monomers are not allowed to enter the colloids. All lengths are measured in the units of the mean radius of gyration of free (ideal) polymers R g ≡ N b 2 /6. We consider a system consisting of n t = 20000 free polymers in a simulation box of size L x = L y = 8, L z = 32, resulting in an invariant degree of polymerization [29] √N =ρ density of the free polymers and R e is the mean end-toend distance of free polymers). The system is discretized in cubic cells of side length 0.25, which are used both for the field-theoretic calculations and the evaluation of local monomer densities. Two colloids of radius R g are placed on the centerline x = y = 0 at fixed distance from each other. They are coated with n b graft polymers with either n b = 37 (low grafting density) or n b = 143 (high graft density). The densities of p-chains are calculated using the particle-to-mesh method [24].
To determine a suitable tuning function ∆µ(r), we must first choose a pair of values ∆µ f and ∆µ p , for which a homogeneous bulk system is occupied almost exclusively by f-chains (fields) or p-chains (particles), respectively. A good choice in our system is ∆µ f = −4 and ∆µ p = 1.2. The function ∆µ(r) then interpolates between ∆µ p close to the colloids and ∆µ f far from the colloids. Specifically, we used a step profile, ∆µ(r) = ∆µ p +(∆µ f −∆µ p )Θ(r −r c ) with the Heaviside step function Θ, where r is the distance to the closest colloid and the shell thickness was chosen r c = 2.5R g . Figure 1 shows the density profiles of p-and f-chains along the line x = 0, y = 0, along with a snapshot of the particle chains in the system. One can see that particle chains dominate close to the colloid, while in the bulk region far from the colloid, the polymers are mostly represented by fields. The total volume of the particle region is roughly ∼ 130R 3 g . The system was studied using a Monte Carlo simulation method which includes three types of updating FIG. 2. Effective force fe between colloids, given by the sum of spring force fs (inset) and contact force fc (inset) in unit of k B T Rg as a function of the distance between the two colloids d with n b = 37 (a) and n b = 143 (b), calculated with the hybrid model (lines), and the corresponding pure particlebased model (symbols). The error bars for the hybrid model are comparable to those for the particle model. steps: (I) update the particle configurations, (II) update the fields using a dynamic density functional scheme, (III) update the {τ } configurations and switch the chain identities accordingly. Moves (I) and (III) are accepted according to the appropriate Metropolis criterion [40]. To assess the performance of the hybrid model, we have also carried out reference simulations of the same system in pure particle representation. The hybrid simulations were roughly three times faster than the simulations of the particle model.
We first consider the effective force between colloid particles [41][42][43], which determines the stability and uniformity [44] of the composite material. It is given by the mean total force acting on one colloid if the other one is kept fixed at a certain distance, and it has two contributions: The mean spring force from the graft polymers, and the mean contact force due to the unsymmetrical collisions of the beads around the colloid. The latter can be expressed [45,46] as an integral over the surface A of the colloid f c = − d 2 A n ρ(A), where n is the surface normal, and ρ(A) the local density of beads at the surface, which includes p-, f-, and graft chains. Figure  2 shows the effective total force as a function of the distance d between the two colloids for different numbers of brush polymers. At low grafting density, the colloids attract each other due to the depletion effect. At high grafting density, the brushes induce an entropic repulsion. This is the regime where the brush stabilizes the colloidal system. For comparison, we also show the results for the reference pure particle system. They are in good agreement with the results from the hybrid model.
Next we investigate how the colloids perturb the surrounding polymer medium. Since the χ parameter (χN = 9) is below the order-disorder transition (ODT) point ((χN ) ODT 10.5 [47,48]), the polymer melt is homogeneous in the bulk. Close to the colloid surface, we observe colloid induced ordering. Figure 3 shows the density profiles for all A-beads, all B-beads and the total density along the line x = 0, y = 0 for a systems containing two at positions (0,0,-11) and (0,0,11), respectively. Only the density profile in half the system is shown, since the other half is symmetric. The results obtained from the pure particle model, also shown in Fig. 3, are again in good agreement, except for a small density dip in the p-f interfacial region. A similar density dip, with comparable magnitude, has also been found in other adaptive resolution schemes [9,10]. In our case, it can be related to the mean-field approximation: When increasing the density of free polymers, the dip becomes smaller (see inset in Fig. 3). It can be reduced by making the "interfacial region" between pand f-regions broader, e.g., choosing a smooth tanh-like profile for ∆µ(r) instead of the simple step function used here. A detailed analysis of these effects will be published elsewhere. Another possibility is to follow Ref. [10] and introduce an additional potential in the interfacial region.
By using a sharply varying tuning function that produces a relatively pronounced density dip, we can assess its influence on the other structural properties of interest. Despite the artifact, colloidal forces are still reproduced accurately by the hybrid model, and the relative distribution of A and B monomers around the colloid is in good agreement with that in the pure particle model. Thus the presence of the artifact seems acceptable in the present system. It might cause problems if one adds small molecules, which might accumulate at the p-f "interface" and whose transport properties across the interface might be altered. In such simulations, the artifact should be removed, e.g., by choosing a tuning function that varies sufficiently slowly.

IV. SUMMARY
In summary, we have developed a hybrid particle-field scheme for simulations with adaptive resolution, which dynamically couples finer particle degrees of freedom with coarser field degrees of freedom. The scheme has been tested at the example of a nanocolloid-polymer composite and verified by comparing results from hybrid simulations to results from pure particle simulations. The new scheme has been derived using a field-theoretic methodology that can be applied very generally to molecular systems without hard core interactions. Hence the approach should be widely applicable for all materials which can be described by Hamiltonians with soft interactions, i.e., typically soft matter systems.
In the present application, the hybrid simulations were found to be roughly three times faster than the corresponding pure particle simulations. The speedup will be even bigger in large systems containing only small regions where a particle representation is necessary. Field-based simulations have the advantage that the computational costs do not increase with the number of molecules. The hybrid approach will thus be particularly attractive for simulations of dense systems, or of polymers with large polymerization index, where particle simulations become expensive compared to field-based simulations. Compared to pure field-based simulations, the hybrid simulation method has the advantage that inclusions and surfaces can be modeled accurately without having to resort to approximate effective descriptions [13].
Since we have focused on equilibrium static properties in this work, we have used a Monte Carlo simulation method to sample the partition function. More realistic dynamical models can be implemented as well. For example, overdamped Brownian particle simulations can be combined in a straightforward manner with a dynamic density functional that reproduces Rouse dynamics in field-based simulations [35]. This model would however neglect hydrodynamic interactions. In order to include these, one could combine a molecular dynamics scheme for the particles [17] with a momentum-conserving fieldbased simulation scheme [49,50]. Such an approach would allow one to use the hybrid model for studying dynamics and flow phenomena in complex fluids.
Another promising direction for future developments will be to replace the tuning function ∆µ(r) that controls the local particle and field content by a function that depends on local densities or order parameters, ∆µ(ρ(r)). The high resolution regimes can then adjust on the fly to the local configurations.
Appendix A: Construction of the hybrid particle-continuum scheme: Technical details For simplicity, we derive the hybrid scheme for a simple polymeric system of n Gaussian chains (labeled α) of one (chemical) type with N monomers (labeled j) in a volume V at temperature T . In the following, energies are given in units of 1/k B T , and lengths in units of the radius of gyration of ideal chains, R g = N b 2 /6, where b is the statistical segment length. The total energy H includes the Gaussian spring energy where R αj denotes the position of monomer j in chain α, and non-bonded contributions described by an Edwards term that is defined in terms of local densities, e.g., with excluded volume parameter v > 0. Here ρ 0φ = α,j δ(r − R α,j ) is the configuration dependent density, which includes contributions from all monomers j of chains α at positions R α,j , and ρ 0 = nN/V is the mean density. The total partition function is then given by In the first step, we partition all chains into two different species, named p-chains and f-chains. This is done by attaching an additional variable (label) τ αj ∈ {0, 1} to each monomer. This spin like variable τ can be coupled to the tuning function ∆µ(r) by exploiting the identity 1 τ =0 exp τ ∆µ(r) − ln e ∆µ(r) + 1 = 1. (A4) This identity holds for any form of ∆µ(r) at any position r, so our method is not restricted to some specific forms of ∆µ(r). Inserting this identity for each τ α,j into the partition function, Eq. (A3), one gets with (A6) where we have defined U ∆µ (r) := ln(e ∆µ(r) +1). This partition function describes a system with additional auxiliary degrees of freedom τ α,j , which however have no physical meaning. The construction ensures that the physics is not changed, compared to the original system.
Let us now assume that we have a given partitioning of molecules into two virtual identities, namely n p p-chains and n f f-chains (with n = n p +n f ). Then the non-bonded energy is given by withφ p the configuration dependent monomer density for the p-chains, andφ f analogously for the f-chains. Obviously, the value of H nb for a given system configuration does not depend on the partitioning into p-and f-chains. Therefore, one can use a different partitioning for each set of {τ α,j }. We shall use the rule that a chain α is an f-chain if τ α,j = 0 for all monomers j, otherwise it is a p-chain. The function ∆µ(r) then tunes the statistical weight of a particular p-and f-chain partitioning in the partition function. The p-chains and f-chains act like two different species, thus the system becomes semigrand-canonical. The procedure so far has only complicated the notation. However, the p-chains and the f-chains can now be treated by different representations. We keep the particle description for the p-chains, but convert the description of f-chains into a field representation. This is done technically in the usual way by inserting an identity operator [27] for the local densitiesφ f of the f-chains, where φ f now denotes the associated density field and iω f is the conjugate field. The particle degrees of freedom of the f-chains can be integrated out resulting in a single chain partition and their associated physics is described by the fields φ f and ω f . In Eq. (A9), a normalization factor has been included for numerical convenience. The partition function then can be written in the form where the index α p indicates that the configurational integral {αp,j} dR αpj now runs over the monomers of p-chains only. The effective Hamiltonian H eff for a given configuration {τ αp,j } with particle positions {R αp,j } and field values φ f and ω f can be split into three contributions Here H p corresponds to the pure contributions of pchains, including in particular their interaction with the virtual potentials ∆µ(r) and U ∆µ (r), H f describes the pure contribution of f-chains in field representation, and finally, the coupling term is given by The partition function given by Eq. (A11) is our final, and formally exact expression of the partition function for the present hybrid particle-continuum scheme. This partition function contains both the particle and continuous field degrees of freedom. Unfortunately, the partition (A11) cannot be sampled efficiently due to the imaginary contribution of iω f , which creates a sign problem (an oscillating integrand). This problem is well-known in field-theoretic polymer simulations [31]. It can be overcome by using the (computationally expensive) Complex Langevin (CL) simulation method [26], but this comes at the expense of having to introduce complex density fields. Hence combining the CL method with particle simulations is not straightforward.
However, most field-based simulation methods operate with real density fields, which is made possible by employing additional (mean-field) approximations. For example, in binary polymer blends, the main effect of fluctuations was found to be sampled correctly by an approach which treats the integral over ω f fields by a saddle point integral and just samples the densities φ f [32]. This can be done within a suitable dynamic density functional scheme [29]. In the present work, we go one step further and also neglect the fluctuations of φ f by setting the noise in the dynamic density functional equations to zero. Such a treatment is known to become accurate in the limit of high polymer densities, or high invariant degree of polymerization [32]. Our approach should be efficient in simulations where large parts of the simulation volume can be treated safely at the (dynamic) mean field level. For example, in phase separated polymer solutions, regions with high polymer densities can be treated at the field level, and regions with low densities at the particle level. chains were generated randomly with Gaussian distributed bonds. More sophisticated schemes such as configurational bias Monte Carlo moves [40] are conceivable as well. Trial moves are accepted or rejected according to a Metropolis criterion. Note that the field ω f and thus the propagator Q f remain fixed in this step.
In our simulations, one "Monte Carlo step" included on average one trial move of R αj per (particle) monomer ((α, j), 2000 trial switches of a variable τ αj (corresponding to one attempted p-f switch per ten chains in our system of 20, 000 chains), and the fields were updated every third Monte Carlo step.