Development of a Meshless Kernel-Based Scheme for Particle-Field Brownian Dynamics Simulations

We develop a meshless discretization scheme for particle-field Brownian dynamics simulations. The density is assigned on the particle level using a weighting kernel with finite support. The system’s free energy density is derived from an equation of state (EoS) and includes a square gradient term. The numerical stability of the scheme is evaluated in terms of reproducing the thermodynamics (equilibrium density and compressibility) and dynamics (diffusion coefficient) of homogeneous samples. Using a reduced description to simplify our analysis, we find that numerical stability depends strictly on reduced reference compressibility, kernel range, time step in relation to the friction factor, and reduced external pressure, the latter being relevant under isobaric conditions. Appropriate parametrization yields precise thermodynamics, further improved through a simple renormalization protocol. The dynamics can be restored exactly through a trivial manipulation of the time step and friction coefficient. A semiempirical formula for the upper bound on the time step is derived, which takes into account variations in compressibility, friction factor, and kernel range. We test the scheme on realistic mesoscopic models of fluids, involving both simple (Helfand) and more sophisticated (Sanchez–Lacombe) equations of state.

In particle-field simulations, the bonded interactions are often described by effective potential energy functions (analogous to those invoked in conventional particle-based schemes), whereas the nonbonded interactions are described by a free energy functional involving a spatial integral of a free energy density depending on one or more particle density fields.The evaluation of the nonbonded free energy entails discretizing the simulation domain to account for the local densities and perform spatial integration.
A popular approach for discretizing the simulation domain is the imposition of a mesh across it.The local density is described at the level of the mesh cells based on the contributions of the particles participating in them. 29,33−44 Because the density is conserved at the cell and not at the particle level, mesh-based approaches can become cumbersome in some cases.For example, when applied to systems with spherical geometry such as droplets and spherical cavities, they may introduce discretization artifacts 40,41 and induce a limit to the maximum resolution.It should be mentioned, however, that recent years have seen the rise of advanced reciprocal-space approaches that provide superior control over discretization artifacts and resolution limitations. 45,46n alternative class of discretization schemes is the so-called meshless approaches.Unlike mesh-based methods that rely on a fixed grid, meshless methods define the domain using the particles themselves.−53 While grid-based methods excel in homogeneous systems (and in many cases perform better), meshless methods excel in simulating large deformations and moving boundaries, e.g., during fracture phenomena. 8,9,47−57 The density field at each particle is estimated by imposing a weighting kernel, which accounts for the contribution of the particles in the local vicinity.It is well known that the kernels invoked by these schemes suffer from various deficiencies, such as tension instability and artificial clumping, 8,47,57−59 and boundary issues. 47Significant effort has been made to address the aforementioned issues in regard to improved weighting kernels, 50,60−62 the so-called artificial stress and viscosity, 9,63,64 staggered meshes and stress-points (stressparticle SPH), 52,65 the moving least-squares SPH, 58,66 the particle shifting scheme, 67 square gradient terms, 13,44,68 invoking ghost particles, 50,69,70 and modifying the weighting kernel 27,44,50−53 at the boundaries.The aforementioned remedies are not crucial in homogeneous systems, as tension instability increases significantly as the system moves far away from equilibrium, in the postfracture regime. 64,71ere we develop a meshless discretization scheme for particle-field Brownian dynamics simulations.By following the footsteps of relevant SPH, 26 SDPD, 68 and MDPD 13−15 implementations, the domain discretization is realized by ascribing an effective number density at the position of each particle as a weighted average of mass contributions from the neighboring particles in the close vicinity.Following the discussion on mesh-based and meshless methods, our approach naturally benefits from the advantages we outlined (such as describing complex geometries and large deformations), but may also experience the drawbacks mentioned.To the best of the authors' knowledge, there are limited (if any) mesoscale simulation frameworks incorporating kernel-based discretization schemes alongside Langevin dynamics in the high friction limit.−76 The central focus of the article is to assess the numerical stability (NS) of the meshless scheme in terms of reproducing the proper thermodynamics and dynamics of homogeneous samples.The weighting kernel is described with Lucy's function. 26,54For the purpose of testing the scheme, our primary analysis is conducted by utilizing Helfand's (HFD) 77 equation of state (EoS), for which the equilibrium density and compressibility can be determined exactly from closed-form expressions.−81 Note, however, that HFD EoS can approximate any EoS under bulk conditions around a reference density; hence, our analysis is directly transferable to more elaborate EoS for continuous phases.The latter is corroborated through additional investigations on realistic fluids by employing the Sanchez−Lacombe EoS. 82,83The performance of the model in inhomogeneous geometries and the effect of the aforementioned remedies on thermodynamics will be explored in a subsequent publication.
NS depends on the choice of model parameters such as monomer mass, coarse-graining degree, thermodynamic conditions, monomeric friction coefficient, time step, and the properties of the weighting kernel.Improper parameter choices can have detrimental effects on the thermodynamic and dynamical properties of the samples and the effect might not be apparent on first inspection.For example, specific parameter combinations may reproduce the correct density but overestimate the compressibility by orders of magnitude.
We invoke a reduced description that takes into account the interrelations of the aforementioned parameters and simplifies our analysis considerably.In particular, we find that NS depends strictly on the reduced reference compressibility of EoS, the ideal number of interacting particles within the range of the kernel, the time step in relation to the friction coefficient, and the reduced external pressure, the latter being relevant when enforcing isobaric conditions.We demonstrate that appropriate parametrization allows for precise restoration of thermodynamics.In practical cases, the accuracy can be further enhanced through a simple renormalization process that involves carrying out an extra simulation.The scheme yields precise dynamics for reasonable parameter choices as well, which can be restored in an exact fashion through a trivial manipulation of the time step and friction coefficient.
We provide a semiempirical relation regarding the upper bound for acceptable time steps, above which the scheme becomes numerically unstable.The relation accounts for the effect of varying the EoS compressibility and friction factor exactly, regardless of the choice of the weighting kernel; the effect of the latter is described in an empirical manner.
The impact of the coarse-graining degree on entropy introduces unexpected side effects in regard to the equilibrium density and compressibility.These effects may lead to unintended complications in some setups.We illustrate the possibility to formulate EoS coefficients by taking account of the coarse-graining degree, ensuring an exact replication of the target density and compressibility.This approach can be extended to more elaborate EoSs.
The article is structured as follows: Section 2.1 provides information regarding the reduced description of the model, Section 2.2 reports the mathematical formulation of the model, Section 2.3 discusses the core assumptions underlying kernels, and Section 2.4 conducts a scaling analysis of the equations of motion, which quantifies the effect of reduced compressibility and friction coefficient on the time step.Section 3.1 displays the capability of the scheme to achieve proper thermodynamics and dynamics across the parameter range considered here, Section 3.2 demonstrates the universal behavior in the reduced description and a renormalization scheme for enhanced performance, and Section 3.3 provides a recipe for canceling the effect of coarse-graining degree on equilibrium thermodynamics.Section 3.4 applies the model to realistic fluids composed of coarse-grained particles.Section 4 concludes the manuscript.The appendices report the density and kernel gradients and demonstrate the equivalence of the derived formula for force with common expressions from the literature.

The Journal of Physical Chemistry B
The Supporting Information, which includes refs 39,72,73,84−89, discusses the implementation of the scheme and the contributions of self-and interparticle interactions to force and stress, provides detailed information regarding the simulation protocol, and presents benchmarks for validating the consistency of the reduced description.where N m is the number of monomers (coarse-graining degree) represented by one particle and m m is the monomer mass, e.g., see Figure 1.The corresponding average particle density measured over the whole system is

Model System and Reduced
The thermodynamics is described by an excess Helmholtz free energy of the form which is modeled in terms of a functional integral of the excess free energy density (a ex ) in conjunction with a square gradient term 13,[42][43][44]68 over the system domain ( ), with n = n(r) being the local number density.The particles execute position Langevin dynamics in the high friction limit (Brownian motion) subject to the stochastic equation of motion with F being a conservative force, ξ a stochastic force satisfying the fluctuation−dissipation theorem, and ζ the friction factor for a particle: where ζ m is the monomeric friction coefficient.
Let n r = ρ r /m be a reference particle density and κ r a reference isothermal compressibility.The model parameters will be nondimensionalized in terms of the reference quantities: where m r is defined as the particle mass, σ r is a characteristic length scale dictated by a reference number density n r , and ε r is the thermal energy.Henceforth, the reduced quantities will be indicated by a tilde above them.As a consequence of Buckingham's π-theorem of dimensional analysis, 90 the reduced description entails that ñ= n/n r = ρ/ρ r = ρ̃and nr = m̃= T ̃= 1; hence, the effect of varying κ r , ρ r , m m , N m , and T is lumped to the reduced reference isothermal compressibility: and the reduced friction coefficient The former can be envisioned as the reference compressibility of the sample relative to the compressibility of an ideal gas with particle density n r at temperature T.
2.2.Evaluation of the Excess Helmholtz Free Energy.The discretization of the reduced excess Helmholtz free energy, can be realized by ascribing an effective number density on the particle level as a weighted average of the number of neighboring particles: 12,25,26,57,62 = < n m m w r 1 ( )

The Journal of Physical Chemistry B
where rj k = |rj k |, rj k = rk − rj, m̃k is the mass of the kth particle, and w̃is a weighting kernel with a finite support at rc.It is worth mentioning Trofimov et al.'s 14 weighting scheme that excludes self-interactions (i.e., j ≠ k in eq 13) and recovers the exact density for homogeneous ideal gas fluids.Subsequently, eq 12 becomes where aẽ x,k �aẽ x (nk), and V ̃k�1/nk=m̃k/ρk is the volume ascribed to each particle.The discretization of the free energy functional into explicit particle-based contributions aligns with derivations found in relevant DPD implementations. 13,14For example, compare eq 14 with eqs 3 and 4 in ref 13, and with eq 16 in ref 14.Note that, even though the "thermodynamic" system volume based on the summation of the particle volumes is not guaranteed to sum to V sim , in dense enough systems the two volumes are expected to be similar. 26he contribution of the excess free energy density to interparticle forces is The right-hand side is expressed in terms of the particle excess stress (σẽ x,k �dA ̃ex.k /dV ̃k, A ̃ex.k = V ̃kaẽ x , k ) and particle density derivative (∇ rd i nk = −nk 2 ∇ rd i V ̃k).As we demonstrate in Appendix D, eq 16 is equivalent to the symmetric form reported in the literature. 12,25,47,57,62he contribution of the SG term to interparticle forces is By expressing ∇ rd i V ̃k with respect to the density gradient and substituting ∇ rd i (∇ rd k nk) 2 from eq 66 in Appendix B we get with C ̃jk being a tensor defined in eq 62, Appendix B. For additional information regarding the density gradients and the kernels, the reader is referred to Appendices A−C.
Owing to the pairwise nature of the interactions, the stress tensor can be determined directly from the forces with the Virial formula. 91The implementation of the scheme and the contributions of the self-and interparticle interactions to force and Virial are reported in the Supporting Information S1.
The excess free energy density will be described with the HFD expression: in terms of the reduced reference compressibility and particle density.The expressions for the corresponding pressure and isothermal compressibility are the following: where the term ñon the right-hand side in eqs 20 and 21 constitutes the ideal gas contribution.It is notable that κr equals the reduced excess compressibility at the reference density (κr = κẽ x (nr)).
The equilibrium density nẽ os at a prescribed pressure can be calculated analytically by solving eq 20 for n: and the corresponding isothermal compressibility κẽ os by inputting nẽ os to eq 21.

Underlying Assumptions and Kernel
Renormalization.Let g(r) denote the radial distribution function of a sample.The mean number of particles at distance r from a reference particle equals, N(r) = n sim g(r).The effective number density of the kth particle can be estimated using a mean-field approach based on the following equation: Conventionally, the weighting kernel is normalized as follows: 13,14,26,68 = dr r w r 4 ( ) 1 By treating the fluids as ideal gases (infinite compressibility, = g r lim ( ) 1 r ), by omitting the self-interactions (w(0) → 0), 14 and by taking account of the normalization in eq 24, eq 23 becomes exact: Under these ideal gas conditions, the number of interacting particles within the range of the kernel can be determined according to eq 26.
For nonideal gas fluids, the equivalence becomes less accurate.This limitation arises from the ad hoc normalization in eq 24, which neglects fluid structure since g(r) ≠ 1.As discussed in Trofimov et al., 14 g(r) features a "correlation hole," i.e., a region of depleted density caused by the repulsive forces between particles.This effect is visualized by the g(r) in Figure 2, which reveals the nonuniform structure of nonideal gas fluids and the depletion region that becomes more pronounced with increasing coarse-graining degree.In addition, eq 24 does not The Journal of Physical Chemistry B take into account the self-interactions w(r = 0), which partially compensate for the correlation hole in g(r).These assumptions improve with increasing r c , where bulk contributions dominate eq 23 (g(r) ≈ 1, number of interactions scaling as ∼ r 3 ) and the self-interactions become less important.
The discrepancy can be further suppressed by renormalizing the weighting kernel with a compensating factor C w , accounting for the self-interactions and fluid structure: In principle, C w can be estimated self-consistently in a meanfield manner by applying eqs 25 and 27 to eq 23 and solving for C w : with n sim,Cw being the mean number density from a simulation carried out with a prescribed value of C w .Physically, C w scales the thermodynamic volume 26 (eq 15) to match the system's actual volume.Interestingly, Trofimov et al. 14 use a similar strategy as in eq 27 to restore the volumetric properties of the fluids, in terms of fitting the density kernel with a linear function accounting for the structure of the fluid. 14,15.4.Reduced Units and Numerical Stability (NS).NS of the nonbonding scheme in terms of properly describing the bulk thermodynamics and kinetics depends on several parameters such as the monomer mass (m m ) and the coarse-graining degree (N m ), the parameters of the EOS (e.g., ρ r and κ r for HFD), the thermodynamic conditions (T, P), the time step in relation to the friction factor, and the properties of the kernel function.
As a previous grid-based study 41 reported, finer discretization (lower N m ) leads to a slight increase in residual stress in bulk samples.The increase was minimal (0−4 atm) for κ SG corresponding to realistic surface tension values.Since the SG term has little impact on bulk properties, we will ignore the effect of κ SG , i.e., set κ SG = 0 from now on.
The thermodynamics in the reduced description is independent of particle mass, temperature, and density (m̃= T ̃= ñr = 1), and can be described by κ̃r.The range of the kernel constitutes a numerical parameter and can be expressed in terms of the ideal number of interacting particles at the reference number density (eq 26): The thermodynamics is independent of the friction factor, the latter being a kinetic property.
We can infer the effect of model parameters on the critical time step Δt ̃crit �below which the scheme is numerically stable�by conducting a scaling analysis on the stochastic differential eq (eq 4) whose reduced form can be discretized (e.g., with Euler's method) as follows: The right-hand side involves a stochastic displacement with variance 2Δt /ζ ̃i, with the triplet of random variables . By applying the expression for the excess part of the force (eq 16) and the excess particle stress from the EoS (σẽ x =−P ̃ex , from eq 20), eq 30 can be expanded as follows: We notice that the conservative part of the displacement scales proportionally with 1/κ̃r and Δt /ζ ̃i.In addition, as long as Δt /ζ ̃i is constant, the dynamics is unaffected.The part within the summation depends strictly on the properties of the kernel (i.e., N c or r̃c).For fixed N c , we get a scaling: t r crit (32)   or, in real units, We note that Δt crit is T-independent and scales sublinearly with the coarse-graining degree (∼N m 2/3 ).

Evaluating the Numerical Stability across a
Reduced Parameter Space.The calculations were conducted in the isothermal−isobaric ensemble (NPT) by following the simulation protocol reported in Supporting Information S2.Throughout the simulations, ζ m , m m , ρ r , N, P, and T were fixed and N c , N m , and κ r were varied (see Table 1).Note the   Interestingly, maintaining the system pressure at a nonzero value (e.g., NPT simulations) introduces an unexpected side effect in terms of NS, in that the external pressure P ̃in the reduced description increases proportionally with N m , i.e., P ̃= Pσ r 3 /ε r ∼ N m .It will be shown that varying P ̃does have an effect on the NS of the scheme, albeit it is relatively weak in our case.It can, however, become considerable with increasing P.
We will explore the NS of the scheme over a broad parameter space in terms of varying N m = {32, 128, 512}, N c = {32, 64 and 128}, and κr m = {0.1, 1, 10}; the latter is the reference compressibility of the sample relative to the compressibility of an ideal gas with monomer density n m at temperature T: NS will be quantified in terms of achieving the equilibrium particle density (ns im ) and isothermal compressibility (κs im ), with respect to their exact values nẽ os (eq 22) and κẽ os (eq 21), respectively.The effect on kinetics will be assessed in terms of comparing the reduced self-diffusion coefficient (D ̃sim ) relative to the exact value from Einstein's model: 93  Increasing N m , κr m , and N c allows working with larger time steps.By taking into account the scaling of N m and κr from eq 32 and optimizing for N c , we derive the following semiempirical formula:  (37)   in real units.In conducting this fit, we have disregarded parameter combinations that result in very poor performance even for small time steps; all cases where N c = 32, and the case where N c = 64, N m = 32, and κr m = 10.Equation 36 yields a reasonable high limit for acceptable time steps as illustrated by the evaluations (vertical dashed lines) in Figure 3. Figure 4 illustrates numerical evaluations of the simulated isothermal compressibility that was estimated via finite differences by conducting an additional simulation at a slightly higher pressure: where P ̃1ε r /σ r 3 = 1 atm and P ̃2ε r /σ r 3 = 5 atm.Similar to density, the compressibility plateaus (diverges) at short (long) time steps commensurate with Δt ̃crit (vertical lines) from eq 36.In terms of reproducing the correct compressibility, the performance is improved with increasing N m (left to right) but most importantly with increasing N c .In particular, setting N c to 32 (64) yields a compressibility, which is about 1−2 orders of magnitude (2 times) higher than κẽ os .On the contrary, for N c = 128, the compressibility is reproduced quite accurately.
The effect of the model parameters in terms of achieving the proper dynamics is illustrated in Figure 5, which presents evaluations of D ̃sim /D ̃Einstein across the parameter space.For low N c , the diffusion coefficient is significantly lower than the predicted value from Einstein's model due to perturbations from  The Journal of Physical Chemistry B the intermolecular interactions.With increasing N c , however, the two become very similar, i.e., see the last row of Figure 5.
Excluding cases that result in very low densities, we notice that the correspondence is improved with increasing compressibility because the sample softens and the perturbations become weaker.The properties of each sample depend strictly on N c and κr, and not on the individual values of κr m and N m .However, because the simulations are realized in the NPT ensemble with P = 1 atm, we have to take into account the effect of the N mdependent external reduced pressure: P ̃= Pσ r 3 /ε r ∼ N m .Beginning our inspection with the most numerically stable case (N c = 128, green markers in Figure 6a−c), we note that the markers collapse on a single curve.The ñs im /ñe os and κ̃s im −1 / κ̃e os −1 exhibit similar behavior: they plateau at ∼1 at low κ̃r and drop monotonically at high κ̃r.On the contrary, D ̃sim /D ̃Einstein increases with κ̃r because the sample becomes softer and the Brownian motion becomes less perturbed.

Universal Response and Renormalization. The left panels of
With decreasing N c , however, there are noticeable deviations from the exact behavior.For example, for N c = 32, the (larger) points that correspond to the largest N m considered here feature increased ns im , κs im −1 , and D ̃sim relative to the other points.This happens because, for low N c , the sample becomes very compressible and as a result the change to the external pressure with increasing N m has a significant effect.
It becomes apparent that the assumptions invoked in the definition of the weighting kernels (see Section 2.3) become very approximate in situations where the range of the kernel is short.Nevertheless, by reweighting the kernel, it is possible to enhance its performance in terms of achieving the correct density, or conversely, to bring the "thermodynamic" volume (V therm , eq 15) closer to the system volume V sim .
The right-hand side panels of Figure 6d−f illustrate the same master plots but from simulations with renormalized weighting kernels, where the C w factor in eq 27 has been optimized with a Newton−Raphson method 94 to reproduce the correct density.The optimized coefficients in each case are displayed in Figure 6a with crosses.By construction, the simulations with the optimized coefficients yield very accurate density, i.e., ns im,Cw,opt ∼ nẽ os in Figure 6d.
From a practical viewpoint, it is notable that in cases where N c and κ̃− 1 are not too low, the optimized coefficients are very close to the ratio ns im /nẽ os from a simulation with the normal kernel (Figure 6a): This is because ns im /nẽ os captures the discrepancy between V therm and V sys , and thus, by multiplying the kernel with ns im /nẽ os , we can partially negate this effect.
The compressibility of the sample with the renormalized kernel improves as well (e.g., compare panels b and e of Figure 6), albeit it remains still much lower than κ̃e os in cases with low and moderate N c .It is thus imperative to set N c to a high value in applications where a correct description of the compressibility is important.
D ̃sim improves slightly in situations where the kernel is renormalized (e.g., compare panels c and f of Figure 6).However, this is of minor importance, since it is trivial to restore the exact (or any arbitrary  In doing so, the diffusivity is reproduced exactly, regardless of whether the kernel is renormalized (see filled markers in Figure 6c,f).

Canceling the Effect of Coarse Graining on
Thermodynamics.Let ρ t and κ t be the target mass density and compressibility of a system with coarse-grained particles, each one having a mass N m m m at temperature T and pressure P. By expressing the free energy density with the HFD model, eq 19, we get the following equations for pressure and compressibility: i k j j j j j y Because the ideal gas contribution (right-hand side of eqs 42 and 43) is a function of the coarse-graining degree ( ), the equilibrium density and compressibility vary with N m .In many applications, this may be an unwanted behavior.
Conveniently, we can calculate analytically the reference density (ρ r = n r N m m m ) and compressibility (κ r ), which reproduce ρ t and κ t as a function of T, P, and N m .First, we solve eq 43 for κ r : By substituting κ r from eq 44 into eq 42, we get the following expression for n r : Finally, we derive the expression for κ r by substituting n r from eq 45 to eq 44: The resulting n r and κ r allow to reproduce the target density and compressibility of the system regardless of the choice of the coarse-graining degree.
3.4.Application to Real Fluids Composed of Coarse-Grained Particles.In this section, we apply the nonbonding scheme for realistic fluids composed of coarse-grained particles, employing various EoS.
Figure 7 presents the equilibrium density and inverse compressibility of coarse-grained water particles at T = 277 K and P = 1 atm, considering different degrees of coarse graining (N m ).The experimental density and compressibility under these conditions are ρ exp = 1000 kg/m 3 and κ exp ∼ 0.5 GPa −1 , 95 respectively.In addition, the self-diffusion coefficient of water molecules is D exp = 1.261 × 10 −9 m 2 /s. 96The calculations are conducted in the NPT ensemble using HFD EoS and renormalized weighting kernels with parameters from Table 2.
The renormalized kernels reproduce the EoS density, as indicated by the filled and textured bars in Figure 7a, which are indistinguishable.On the contrary, the inverse compressibility is slightly lower than its exact value κ eos −1 .The optimal C w for each case is in good match with its estimation from the mean-field formula in eq 28 (see the caption of Figure 2).The effect of C w on g(r) is negligible (compare the colored lines with the thin black lines in Figure 2).In addition, the renormalized friction factor reproduces the experimental diffusion coefficient exactly.
Cases W 1 , W 10 , and W 100 in Figure 7 (where the subscript in this notation denotes N m ) illustrate a naive approach where the EoS coefficients are set to experimental density (ρ r = ρ exp ) and compressibility (κ r = κ exp ).For N m = 1 (system W 1 ), the density is reproduced poorly due to the contribution of the ideal gas term (see discussion in Section 3.3).As N m increases, however, the discrepancy progressively decreases.In cases W 1t , W 10t , and W 100t , the coefficients ρ r and κ r have been determined from eqs 45 and 46, respectively, ensuring equality between EoS density and compressibility with their experimental counterparts.
Figure 8 illustrates the impact of density variation with pressure for systems W 1t and W 100t , comparing approaches with and without kernel renormalization (C W = 1).The renormalized  The Journal of Physical Chemistry B kernels lead to volumetric properties that closely match the analytical predictions from EoS in eq 42.Notably, the bulk density at 1 atm is reproduced perfectly, while the pressure− density slope is in good agreement, albeit slightly weaker due to the slightly higher compressibility (see Figure 7b).Conversely, evaluations using non-normalized kernels exhibit a significant deviation from the analytical predictions, with this discrepancy becoming more pronounced with decreasing N m .Intuitively, increasing N m should hinder the accuracy of the model, as it makes the correlation hole in the g(r) more pronounced (e.g., as shown in Figure 2).However, the selfinteractions become more effective at compensating for the correlation hole with increasing N m , thus leading to improved model behavior and compensating factors closer to unity.The dot-dashed line in Figure 8 schematically depicts the relationship between the density of the sample at P = 1 atm using the non-normalized kernel and the optimal C w for each case; according to eq 38, the optimal C w is close to the ratio of ρ sim (C w = 1)/ρ eos .As shown in Figure 8b, the sample with N m = 100 cannot withstand large negative pressures, indicating that susceptibility to cavitation increases with the coarse-graining degree.
Figure 9 showcases applications in homogeneous polymer melts described by the Sanchez−Lacombe EoS 82,83 with excess pressure:  with parameters the characteristic temperature (T*), pressure (P*), and density (ρ* = mn*).For simplicity, the particle molar mass is set to M = 1000 g/mol, assumed to be equal to the polymer molar mass.Note that the effect of chain length is naturally accounted for by the Sanchez−Lacombe EoS through the ideal gas term.The latter is typically expressed as P ig = P*ρT/ (ρ*T*r), with r = P*M/(ρ*RT*) being the number of lattice sites occupied by the polymer chain.
For each case, the reference density (ρ r = mn r ) is determined by solving eq 47 for n at zero pressure.Subsequently, the reference compressibility is set to κ r = κ ex (n r ).Similar to previous cases, the kernel is renormalized through a single simulation, resulting in a simulated density that aligns perfectly with ρ eos , i.e., the filled and textured bars in Figure 9a are indistinguishable.Moreover, the simulated compressibility is slightly underestimated, albeit the effect can be suppressed by further expanding the range of the kernel.
In many practical scenarios, it is convenient to model the polymer chains in terms of multiple particles connected by entropic springs.To achieve the correct scaling of dynamics with chain size and reproduce the experimental viscoelastic properties, 5,33,98 the scheme should be employed in conjunction with a model that accounts for the effect of entanglements. 34,35Various models, such as TWENTANGLEMENT, 99 slip-link, 36,37 and slip-spring 5,33,38,100 models developed in the literature, address this aspect effectively Table 3.

CONCLUSIONS
The article develops a meshless discretization scheme for particle-field Brownian dynamics simulations.The density is ascribed on the particle level in terms of imposing a weighting kernel with a finite support (r c ).The free energy of the system is described by a free energy density from an EoS in conjunction with a square gradient (SG) term.
The contributions of the EoS and the SG terms to force and stress are derived analytically by differentiating the free energy.We demonstrate that the resulting expression for the EoS contribution to force is equivalent to the common expression used in the literature.
The focus of the article is on determining the numerical stability (NS) of the scheme in bulk conditions.In doing so, we primarily invoke Helfand's (HFD) EoS with parameters the reference density and compressibility, and we drop the SG term.HFD EoS can serve as an approximation for any EOS under bulk conditions around a reference density.Our analysis is directly   3.

The Journal of Physical Chemistry B
applicable to more sophisticated EOS and this has been corroborated here by simulations of realistic fluids utilizing the Sanchez−Lacombe EoS.
The NS of the scheme depends on several parameters, including the monomer mass and coarse-graining degree, the coefficients of the EoS, the thermodynamic conditions (P, T), the time step, the friction factor, and the properties of the weighting kernel.To take into account the interrelations of the aforementioned parameters to the NS, we invoke a reduced description of the model, which simplifies the analysis considerably.
We find that the outcome of the simulations depends strictly on the reduced reference compressibility (κ̃r), the ratio of the friction coefficient to the time step (ζ ̃/Δt ), the range of the kernel function (rc), and the reduced external pressure (P ̃), the latter being relevant in isothermal−isobaric statistical ensembles.
NS is assessed in terms of reproducing the equilibrium density and compressibility dictated by EoS and the self-diffusion coefficient from Einstein's model. 93For large enough rc, the aforementioned properties are reproduced very accurately.With decreasing kernel range, on the other hand, we observe significant discrepancies because the assumptions invoked in the definition of the kernels become inaccurate.
We develop a scheme for improving the aforementioned deficiencies in terms of renormalizing the weighting kernels.The optimal renormalization of the kernel requires conducting simulations in the NPT ensemble, to assess the discrepancy of the density (ns im /nẽ os ) and diffusion coefficient (D ̃sim /D ̃Einstein ) relative to their exact values.In practical cases with sufficiently high rc, the density can be restored by multiplying the kernel with ñs im /ñe os from a single simulation.In addition, given that D ̃Einstein ≡ 1/ζ ̃and that the dynamics depends on ζ ̃/Δt ̃and not on the individual values of ζ ̃and Δt , it is trivial to restore the exact diffusion coefficient just by multiplying ζ ̃and Δt with D ̃sim / D ̃Einstein .Restoring the compressibility exactly is, however, not trivial; therefore, in applications that necessitate high accuracy, the usage of kernels with long enough range is advised.
The parameters of the model have a strong effect on the choice of acceptable time steps.The effect of κr and ζ ̃on the time step was determined exactly by conducting a scaling analysis on the equations of motion.Even though it is not straightforward to assess the effect of rc analytically, the latter effect was fitted to a semiempirical formula, eqs 36 and 37, which yields very reasonable estimates regarding the upper bound of acceptable time steps.
The equilibrium density and compressibility depend on the reference density and compressibility of the EoS, the thermodynamic conditions (P, T), and the coarse-graining degree (N m ).In many applications, the latter introduces undesirable side effects but, as we demonstrate, it is possible to derive analytically N m -dependent reference densities and compressibilities that reproduce the target equilibrium density and compressibility exactly.
As will be demonstrated in a future publication, the discretization scheme is compatible with the slip-spring BD/ kMC model developed by the authors, which has been applied for investigating the rheology of linear 5,33,101 and branched 102 high molar mass polymer chains, the elastic properties of rubber materials, 39,98 as well as interfaces between molten polymers and gases 41 or solids. 40uture directions of this study include the application of the kernel-based discretization scheme to mesoscopic polymer/ solid and polymer/polymer interfaces of entangled, high molar mass polymer chains under quiescent and flow conditions, the assessment of the capability of the scheme to predict equilibrium morphologies and interfacial free energies in conjunction with the square gradient theory, and applications in material fracture.

■ APPENDIX A VECTOR GRADIENTS
Let r jk = r k − r j be a vector pointing from r j to r k .The length of the vector is denoted with = || || r r jk jk , and the unit vector with rĵ k =r jk /r jk .The gradient of the vector length r jk with respect to a position vector r i is = r r ( ) jk jk ki ji r i (49)   with δ ij being the Kronecker delta function; δ ij = 1 if i = j and δ ij = 0 is i ≠ j.The partial derivative of a vector r jk with respect to the component α i ∈ (x,y,z) of a position vector r i is The partial derivative of a component αĵ k of a unit vector rĵ k with respect to the component β i of a position vector r i is The gradient of the unit vector rĵ k with respect to a position vector r i yields the Jacobian:  By evaluating the partial derivatives of the unit vector components with eq 51, we derive the compact form:

The Journal of Physical Chemistry B
with r̂j k ⊗r̂j k being a dyadic tensor and I 3 being the unit secondorder tensor in three dimensions.

Gradient of the Number Density
The gradient of the number density at the kth particle with respect to the position vector of the ith particle (eq 13) is evaluated in conjunction with eq 49 as follows: (54)   For i = k: For i ≠ k: (56)   Note the relations between ∇ rd k n k and ∇ rd j n k : and ∇ rd i n k and ∇ rd k n i , for k ≠ i: Equation 58 is a demonstration of global momentum conservation. 58

Double Mixed Gradient of Number Density
The double mixed gradient of the number density can be determined by taking the gradient of eq 54: By substituting ∇ rd i rĵ k from eq 53 and ∇ rd i r jk from eq 49 we get i k j j j j j j j j j j y { z z z z z z z z z z i k j j j j j j j j j i k j j j j j j j j y  (60)   or in a more compact form:

Gradient of the Square Gradient
The gradient of the square gradient can be determined as follows: The first and second derivatives are the following: The expressions of A jk and B jk for calculating C ij in eq 62 are the following:  which is the expression commonly used in the literature. 12,25,47,57,62Note that the summation excludes the reference particle, 25 even though it might not always be explicitly stated.
Description.Consider a system with N coarse-grained (CG) particles occupying volume V sim at temperature T, each of mass =

Figure 1 .
Figure 1.Schematic illustration of a fluid composed of N coarse-grained particles.Each particle represents a set of N m monomers, each characterized by mass m m and monomeric friction coefficient ζ m .The system is simulated in the NPT ensemble under atmospheric pressure.Particle interactions are governed by an EoS subject to a reference density (ρ r ) and compressibility (κ r ).Particle density (n k = ρ k /m k ) is estimated using a weighting kernel (w) that takes into account the influence of neighboring particles within a specified cutoff distance r c .

Figure 2 .
Figure 2. Radial distribution function g(r) versus r/r c from systems W 1t (solid line, r c = 12.2 Å), W 10t (dashed line, r c = 26.3Å), and W 100t (dotted line, r c = 56.8Å) in Table2.The thin black lines correspond to results for C w = 1.The dashed horizontal line displays the g(r) of an ideal gas.The compensating factor C w from the numerical fitting (from the self-consistent relation in eq 28) equals 0.984 (0.982), 0.994 (0.993), and 0.998 (0.997) for W 1t , W 10t , and W 100t , respectively.

Figure 3
Figure 3 illustrates the equilibrium particle density as a function of N m (left to right), N c (top to bottom), and κ̃r m

Figure 7 .
Figure 7. (a) Density and (b) inverse compressibility of coarse-grained water particles using the parameters listed in Table 2.The colored bars illustrate the simulated quantities, while the textured ones indicate the exact values determined by EoS, i.e., ρ eos and κ eos from Table 2.The dashed lines correspond to the experimental density and inverse compressibility.

Figure 8 .
Figure 8. Pressure as a function of density for systems (a) W 1t and (b) W 100t from Table 2, with ( □ ) and without (○, C w = 1) kernel renormalization.These simulations were conducted in the canonical ensemble (NVT).The black dashed line displays the theoretical prediction based on eq 42.The blue dot-dashed line denotes the densities C w ρ sim,1 atm equal to (a) 984 kg/m 3 and (b) 998 kg/m 3 .The dotted lines are guides to the eye.

Figure 9 .
Figure 9. Same as Figure 7 using the parameters listed in Table3.

Table 1 .
Simulation Parameters in Real and Reduced Units.Note that ζ m , m m , and ρ r Correspond to the Monomeric Friction Coefficient, Monomer Mass, and Density of High-Density Polyethylene (HDPE) at 450 K 5,92 correspondence between mass (m) and molar mass (M), measured in kg/mol: m = M/(1000N A ) with N A being Avogadro's number.As discussed in Section 2.4, arbitrary combinations of Δt, ζ m , m m , ρ r , N m , and κ r which yield the same κr and Δt /ζ ̃i result in the same reduced sample properties, e.g., see benchmarks in the Supporting Information S3.As our analysis focuses on the bulk properties of the samples, the square gradient term (κ SG ) is set to zero in this study.