Comparison of equilibrium techniques for the viscosity calculation from DPD simulations †

Dissipative Particle Dynamics (DPD) is a powerful mesoscopic modelling technique that is routinely used to predict complex fluid morphology and structural properties. While its ability to quickly scan the conformational space is well known, it is unclear if DPD can correctly calculate the viscosity of complex fluids. In this work, we estimate the viscosity of several unentangled polymer solutions using both the Einstein and Green–Kubo formulas. For this purpose, an Einstein relation is derived analogous to the revised Green–Kubo formula suggested by Jung and Schmid, J. Chem. Phys. , 2016, 144 , 204104. We show that the DPD simulations reproduce the dynamical behaviour predicted by the theory irrespectively of the values of the conservative and friction parameters used and estimate a Schmidt number compatible to that of a fluid system. Moreover, we observe that the Einstein method requires shorter trajectories to achieve the same statistical accuracy as the Green–Kubo formula. This work shows that DPD can confidently be used to calculate the viscosity of complex fluids and that the statistical accuracy of short trajectories can be improved by using our revised Einstein formula.


I. Introduction
Dissipative Particle Dynamics (DPD) is a mesoscopic method widely used to predict the morphology of complex fluids in both toy models 1-3 and systematically parametrized models for specific chemical systems. [4][5][6][7] Despite its popularity to quickly scan large shallow energy surfaces, the ability of DPD to predict dynamical properties (e.g. viscosity) is more debated since the methods used to calculate these properties are subject to significant statistical error. 8,9 Non-equilibrium simulations, where shear is applied to the simulation box, are routinely employed in DPD to calculate the system viscosity with minimal computational cost [9][10][11] and have been tested on both simple fluids 8,9,12 and polymers. 13,14 However, it remains unclear whether the equilibrium approaches for the viscosity calculation can be confidently used with DPD. This is a crucial point as one wants to calculate the viscosity of the morphologies resulting from the simulations without perturbing them. Due to the soft inter-bead potential typically employed in DPD models and the high shear rates required by non-equilibrium methods to be computationally tractable, the fluid microstructure is easily lost. 2 Another drawback of the non-equilibrium approaches is that the low shear plateau is not easily achieved in some polymer cases, 12 and the zero-shear viscosity can only be obtained by extrapolating from the lowest shear rate values. 13 Thus, equilibrium techniques are preferred if the zero-shear viscosity is required.
The equilibrium approach uses the integration of the pressure autocorrelation function to calculate viscosity. The major drawback of this approach is the large fluctuations associated with the pressure tensor elements' values, leading to poor accuracy in the resulting viscosity and therefore requiring very long simulations to ensure good statistics in the pressure autocorrelation function. The viscosity resulting from equilibrium simulations is usually calculated using the Green-Kubo (GK) relation. However, the use of the GK relation to calculate bulk viscosity using DPD has been under question since the presence of stochastic force creates extra noise into the resulting stress tensor. 8,14 Moreover, in DPD simulations, uncertainty in the results has also been linked to the use of the dissipative force, where the error in viscosity increases with increasing the friction coefficient. 8 Different integration algorithms 14 and the effect of finite time step on the dissipative and random terms 15 also play a role in the low statistical accuracy. The latter issue, along with the lack of the time reversibility invariance of the dissipative force, were addressed by the recent suggestion of a revised GK formula for DPD simulations proposed by Ernst and Brito 16,17 and implemented by Jung and Schmid. 18 In the revised GK expression, the potential and dissipative forces are integrated, and the random force contribution to the stress is taken outside of GK integral as an additional term.
An additional challenge associated with the equilibrium method is to decide on the limit of the integration point. Jung and Schmid 18 integrate the autocorrelation function up to the point where the curve reaches 1% of its initial value and the remaining values are fitted with a decaying exponential function. In some studies, the entire viscosity integral is fitted by an analytical expression and the infinite analytical limit of this expression is used to obtain the viscosity value. 19 Various analytical expressions for fitting the viscosity integrals have been suggested in these cases, following the work done with conventional molecular dynamic simulations. 20 Both approaches contain a certain level of arbitrariness and lack solid physical reasoning.
An alternative approach for evaluating viscosity is the Einstein formula initially suggested by Helfand, 21 whose equivalence to the GK formula is well known. 22 Later, an intermediate to GK and Einstein formula was suggested to resolve momentum discontinuity for systems with periodic boundaries. 23,24 Furthermore, Español 25 recently presented a generalized Einstein-Helfand form that can be used in the case of underlying stochastic dynamics. Many researchers have adopted this method to calculate the viscosity for Lennard-Jones fluids; [26][27][28][29][30] however, an application of the Einstein viscosity formula in the presence of dissipative and stochastic forces is lacking in the literature to the author's best knowledge. The importance of testing the validity of Einstein-Helfand relations in DPD simulations has also been pointed out in the past, 18 however, this formalism has not yet been applied to DPD systems. One of the main advantages of the Einstein method is that the viscosity can be calculated using the positions and velocities of the particles that are calculated during the numerical integration of the equation of motion. In this way, the unwrapped particles positions and velocities can be used before the periodic boundary conditions are applied, and the momentum discontinuity does not constitute any problems. Moreover, any extra statistical inaccuracy coming from the pressure tensor is avoided. However, for a generalized approach where the viscosity is broken down into the various contributions of the DPD forces, this computational methodology cannot be applied as both positions and velocities are updated using the total force acting on a particle, not the individual force contributions. Therefore, the GK and Einstein formulas' intermediate method is desired in this case.
In this work, we are interested in testing whether DPD simulations can be reliably used to calculate the viscosity of complex fluids from equilibrium simulations. To do this, we simulate various unentangled polymer solutions using a simple toy DPD polymer model with a range of conservative and friction parameters. To calculate the viscosity using the Einstein formalism, we derive an analogous Einstein expression to the revised GK formula suggested by Jung and Schmid. 18 The viscosity values are then compared with scaling laws from theoretical predictions. Finally, we looked at the effect of the conservative and friction coefficients on the viscosity and dynamic behaviour of the polymer systems. We believe that the long-time relaxation of the polymers make them an ideal proxy for structured liquids such as surfactant solutions, for which DPD has previously been shown to be a powerful method to reproduce their phase diagram and whose non-Newtonian viscosity behaviour resembles that of living polymer solutions. 31

II. Methodology
A. Simulation details DPD simulations are run for three different concentrations of polymers c p = {0.3, 0.8, 1.0} and the number of solvent beads in the system is adjusted so that the total number density is r = 3. Four different polymer lengths are studied, N = 10, 25, 40, 50, where N is the number of beads per chain, in a volume of V = L 3 = 20 3 . The number of polymer chains for each concentration is given in Table 1. The FENE (Finitely Extensible Nonlinear Elastic) potential was used for the bond interactions with the literature's parameters. 3 The simulations are performed in the canonical ensemble with target thermal energy k B T = 1.0 where T is the temperature and k B Boltzmann's constant. The effect that the integration algorithm has on properties such as viscosity been investigated in the past 14,32 and those studies were considered here in making the decision of the most appropriate algorithm to use. Since the results indicate that the calculated viscosities fall within range of theoretical predictions and nonequilibrium studies independently on the algorithm used, the standard and most used Groot-Warren Velocity Verlet time integration scheme was employed here 33 to solve the equation of motion with a time step Dt = 0.04. The production runs are 3 million steps long and 10 independent runs are used for the viscosity calculations. We accumulate the pressure tensor every 3Dt. This value was found to give the same results as accumulating every 1Dt whilst reducing the size of the data to be analysed. All of the simulations are performed using the DL_MESO software. 34 Although the equations that describe the DPD forces can be found in the literature, 35 the conservative, dissipative and random forces are given here. The form of the conservative force is with r ij the distance between i,j beads and r c the cut-off radius beyond which no interactions are taking place. The parameter a ij is called the conservative coefficient. The vector e ij is the unit vector defined by the distance between bead i and bead j. Similarly, the dissipative force is defined as with v ij the velocity difference vector v ij = v i Àv j , v i and v j are the velocities of particles i and j respectively and g ij is the friction coefficient.
The stochastic random force compensates for the lost degrees of freedom after coarse-graining and along with the dissipative force acts as a thermostat for the system.
where w ij is a randomly fluctuating variable that follows Gaussian statistics with zero mean and unit variance and ffiffiffiffiffi Dt p the square root of the integration timestep. The random weight function is related to the dissipative weight function through the relation (2.4).
The fluctuation-dissipation theorem relates the g and s parameters with the relation s ij 2 = 2g ij k B T. 36 The DPD beads are topologically connected using the FENE potential with constant , r max = 1.3r c and r o = 0.7r c .
We study six combinations of a ij , g ij . In all cases the conservative and friction coefficients between like-particles are kept constant and are given in Table 2 where the letter P refers to the polymer beads and S to the solvent beads. We consider an a thermal solvent where a ij = 25, g ij = 4.5, two additional a ij (15,20) values at g ij = 4.5 and two additional g ij values (5, 10) at a ij = 25. Varying a ij and g ij corresponds to modifying the solvent quality and bulk solvent viscosity, respectively. The cut-off radius is kept constant at r c = 1 in all the cases and the mass is equal for all beads (m = 1). These two quantities, together with the energy k B T = 1 set the distance, mass, and energy scales respectively.
Since the cut-off radius of the beads (and mass) is the same regardless of the bead type, the polymer concentration in the solution is equivalent to its volume (or mass) fraction, f. Thus, a concentration of c p = 0.3 is the same as volume fraction f = 0.3. From now on, both notations will be used interchangeably to express the polymer concentration and/or volume fraction.

B. Structural and dynamic analysis
We calculated some static properties of the polymer chains, i.e. the end to end distance and radius of gyration, to check the consistency with previously published data. 3 The end-to-end vector, R ee ! , is defined as the vector between the first and last monomer of a polymer chain and its value indicates whether the polymer is in an extended or coiled conformation. 37 The Euclidean norm of ee ! defines the end-to-end distance, R ee . The autocorrelation function of R ee ! gives a measure of the equilibration of the system. Usually, an exponential function of the form of eqn (2.6) is fitted to the autocorrelation function of the R ee ! to calculate the relaxation time of the polymer chain. 38 where A and t r are the fitting coefficients, t r is the relaxation time of R ee ! , and t o the time origin (for an autocorrelation function example, see Fig. S5 in ESI †). We then use t r to determine the cutting point for the viscosity integral. More specifically, the final viscosity values reported below are averaged over the time interval [2t r , 3t r ]. 13 The autocorrelation function in eqn (2.6) is calculated using multiple time origins for accuracy.
Another measure of the expansion or contraction of the polymer chain is the radius of gyration, R g . The R g is calculated as the average over all time frames and monomers of the distance between the chain monomers and centre of mass of the polymer chain. Monitoring R g as a function of the polymer molecular weight (which in the present case corresponds to the number of chain beads, N, since m = 1) provides an indication of the solvent quality.
For the dynamic properties, we calculated the diffusion coefficient, D from the particles mean squared displacement, and the Schmidt number, which is a dimensionless number defined as the ratio between the kinematic viscosity, m, and the diffusion coefficient, 12 In this work, we use the polymer contribution to the solution viscosity, i.e., Z À Z s , for the Schmidt number calculation (eqn (2.8)) where r is the density of the fluid, Z is the viscosity of the solution, Z s is the viscosity of the solvent (Z s = 1.1 in our case, see Fig. S1, ESI †) and D is the diffusion coefficient of the polymer chain. When simulating polymer melts, the viscosity used in eqn (2.8) is the bulk viscosity Z. In this work, D was Table 2 DPD model parameters used in the current study. P and S stand for polymer and solvent beads, respectively. The conservative, a ij , and dissipative,g ij , cross-terms parameters are varied as follows a ij = 15, 20, 25 and g ij = 4.5, 5, 10. In the polymer melts, the parameter of self-particle type interactions is varied as g ii = 4.5, 5, 10 while a ii = 25 always calculated through the mean squared displacement of the chains' centre of mass rather than the displacement of a single chain monomer.

C. Viscosity from Green-Kubo formula
Green-Kubo (GK) relations have been extensively applied to calculate transport coefficients from DPD simulations. 14,16,17,39 The relation that is used for the viscosity calculation is the integral of the pressure autocorrelation function and has the following original form: with P ab being the pressure tensor's off-diagonal ab element with a, b = x, y, z. The symbol h. . .i t o refers to the evaluation of the autocorrelation function over multiple time origins t o .
All the off-diagonal terms of the pressure tensor are used to reduce the statistical error. 38,40 Mondello and Grest have suggested using all nine pressure tensor components by adequately weighting the diagonal elements. 41 Our simulations are extremely long and using only the three off-diagonal elements is sufficient to ensure good statistics.
In DPD, the pressure tensor is the sum of four different contributions: the pressure due to the dissipative force, P D ab , the conservative force, P C ab , the random force, P R ab and the kinetic part of the energy, P K ab . In literature, 42,43 the total pressure resulting from these four contributions is mainly used with the GK relation. However, recently it has been shown that, for high density systems (r 4 4), if all the contributions are included within the integral of eqn (2.9), the resulting viscosity in the linear regime deviates from the value obtained by using non-equilibrium approaches. 18 Thus, a revised GK formula, eqn (2.10), was suggested 16,18 to achieve results comparable to non-equilibrium simulations and resolve the time reversibility issues arising from the dissipative force. 16,17 The random pressure is decorrelating instantly within the first integration time step Dt and thus its integral can be calculated by the basic geometry rule which, with the appro- , gives the instantaneous viscosity Z GK N in eqn (2.10). 18 The instantaneous viscosity term is important for retrieving viscosity values comparable to non-equilibrium methods (see for more details ref. 18). The revised relation (2.10) gives the same viscosity as the non-equilibrium simulations irrespectively to the system density. 18 Therefore, due to the extended applicability of eqn (2.10) over eqn (2.9), the relation (2.10) is used in this work. In the literature 44,45 it has been shown that the highest contribution to the overall pressure autocorrelation function comes from the bonding, the conservative and kinetic forces. Thus, the pressure in eqn 2.10 includes all these terms (i.e. P P ab = P bonds ab + P C ab + P K ab ) and the total viscosity is computed as follows: To estimate the error, we follow the time decomposition method 19 over 10 independent runs. This number of independent trajectories is sufficient to ensure a standard deviation lower than the mean viscosity value. In contrast to ref. 19 the viscosity cut-off point is decided according to the value of the relaxation time of the end-to-end vector, and the average viscosity is averaged over the interval [2t r , 3t r ].

D. Viscosity from Einstein formula
An equivalent expression to the GK relation was derived by Helfand, 21 taking into account the momentum diffusion in the system. It is an Einstein type relation and has the original form: r ia ðtÞp ib ðtÞ where r and p are the position and momentum of the particles, respectively and a,b = x, y, z. To overcome the discontinuity introduced by the periodic boundary conditions, an alternative formula was suggested using the integral of the pressure tensor P ab . 24,46 So, eqn (2.12) becomes: To make a direct comparison with the revised GK formula of eqn (2.10), an equivalent Einstein formula is suggested in this work, which is defined as where each term Z AB is given as in eqn (2.15).
A and B are the potential, P, dissipative, D, or random, R, contributions. If A = B then it is simply the integral squared. 26 The derivation of this formula is straightforward if one starts from the equivalent relation (2.10). For the derivation, the indexes ab indicating the directions are omitted from the general P(t) symbol. However, the derivation holds in all directions. In what follows, the symbol P refers to the total pressure due to the total force contributions. It is known that the GK and Einstein formulas are equivalent and related through the relation: Keeping only the integral terms, I, and substituting them with the Einstein form due to eqn (2.16): The only terms that remain from expression (2.17) are Z PP and Z DD since the terms Z PD and Z DP are equal and opposite and thus eliminate each other. Following Jung's et al. 18 derivation we define the Einstein correction term due to the random force as follows, By multiplying with the volume on each side, we get Alternatively, and in compliance with the Einstein formalism the Z E N term can be also written as in eqn (2.21), Eqn (2.21) together with eqn (2.18) constitute the revised Einstein formalism (eqn (2.14)). As for the GK case, the viscosity was calculated over the time interval [2t r , 3t r ] where the system was well equilibrated, and the statistical accuracy was less than half of the viscosity value. 13

A. Scaling of static properties
In polymer theory, the structural and dynamical properties of polymer chains scale with the number of chain monomers, N, following specific scaling laws. The scaling laws and the scaling coefficient values depend on two factors: the polymer-polymer overlap and solvent quality. 47 For highly dilute solutions, where polymer-polymer overlap can be neglected, the quality of the solvent determines the magnitude of the swelling of the chain as a function of N and the radius of gyration scales with N as where v is the Flory exponent. If the solvent effect is unimportant (i.e., y-conditions) v = 0.5, whilst v E 0.6 indicates that the interactions between the polymer and solvent beads are energetically favourable (i.e., good solvent conditions). Upon increasing the polymer volume fraction, f, above an , Flory screening sets in, and the Flory exponent v again becomes equal to 0.5 for all solvents. Upon increasing polymer concentration, a change in the dynamic properties mirrors the scaling of the structural ones. For dilute solutions, hydrodynamic interactions (HI) are dominant, and the Zimm model describes the scaling with N of the polymer viscosity and diffusion coefficient as Z p N 3nÀ1 and D p N Àn respectively. For highly concentrated solutions and polymer melts, HI are screened, and the dynamical properties of the system follow a Rouse model where Z p N and D p N À1 . 47 The crossover between the dilute and concentrated regimes does not occur suddenly, and for semidilute solutions, the polymer properties scale with the polymer volume fraction. 48 Here we consider the semidilute and concentrated solutions where the polymer volume fraction is above the overlap volume fraction (f*). As the scaling behaviour is fully developed for f f Ã ) 2 3 , we expect to see scaling exponent coefficients in better agreement with the theory in this region although the fitting of the data are reported for 1 o f f Ã o 2 in some of our plots. Jiang and co-workers 3 have shown that for dilute and semidilute athermal polymer solutions, DPD models reproduce the correct scaling of the polymer radius of gyration (R g ) and end-to-end distance (R ee ) as a function of the polymer volume fraction. Here we perform a similar scaling analysis to check the dependence of the Flory parameter value, n, on the DPD interaction parameters. As said above, in semidilute solutions, the scaling of the R g depends on the polymer volume fraction, f, as described by eqn (3.1). 47 where R g0 is the radius of gyration of a polymer chain at infinite dilution. The values of the rescaled R g as a function of the rescaled volume fraction for all conservative and dissipative coefficients at c p = 0.3 and 0.8 are presented in Fig. 1 (left).
Fitting the data with eqn (3.1) we obtain a slope of 0.13 for the two lowest conservative parameters (a sp = 15 and 20), and 0.10 for a sp = 25. Regardless of the conservative interaction parameters, these values are in agreement with the scaling behaviour expected for real chains dissolved in good solvent for which n D 0.59. 47 For a polymer melt, i.e. c p = 1.0, for which excluded volume interactions are thoroughly screened, and the polymer chains follow the Rouse statistics, R g scales as R g E (N À 1) n , where (N À 1) is the number of bonds, with n = 0.5 49,50 ( Fig. 1 (right)). In addition, no evident effect of the g ij parameter on R g is observed. Similar results are obtained calculating the R ee (see ESI † document and Fig. S3, ESI †).

B. Numerical verification of revised Einstein formalism
The two equilibrium viscosity methods are initially tested on a simple DPD ''water'' model comprising of monomers interacting through a ij = 25 and g ij = 4.5 The viscosity for this model was found to be 1.1 using both methods when integrating the pressure tensor to t = 0.8 (see ESI † document and Fig.  S1, ESI †). This result agrees with non-equilibrium DPD simulations carried out using Lees-Edwards boundary conditions, which calculated the zero-shear viscosity to be 1.08. Next, we proceed with the polymer solutions, where ten independent runs are employed to evaluate the error. The average value of the viscosity obtained with both methods is comparable as these methods are numerically equivalent. 22 However, their standard deviation shows a small but significant difference as the polymer length (and relaxation time) increases, with the Einstein method having a clear advantage over the GK one. To compare the two methods, the relative Mean Absolute Error (MAE) of the mean viscosity values given in eqn (3.2) is used.  The symbol h. . .i indicates average taken over all the polymer concentrations and interaction parameters. On average, the viscosity obtained from the Einstein relation had a deviation of B2.4% from the value obtained by GK (eqn (2.9)) ( Fig. 2 and  3). This difference can be attributed to numerical precision errors and does not depend on the interaction parameters. Although the MAE values calculated for both methods show are similar (see Fig. S2 in ESI †), the standard deviations are different. The results indicate that the Einstein method had a lower standard deviation, especially for long polymer chain lengths. Similarly to the MAE for the mean value of Z shown above, we calculate the difference in standard deviations using eqn (3.3).
The absolute value is omitted from the numerator to check whether the standard deviation is smaller in GK (negative MES) or Einstein relation (positive MES). As shown in Fig. 4, for N = 40 and 50, the standard deviation calculated using the Einstein relation can be, on average, 5-8% smaller than the standard deviation obtained using the GK formula. The difference in standard deviations is insignificant for chains with N = 10, 25.
We ascribe this result to the fact that the longer the polymer chain is, the slower is its dynamics and therefore the longer the DPD trajectory needs to be to reduce the statistical error. We observe again that the value of the MES does not depend on the interaction parameters or the polymer concentrations. Our findings are consistent with the literature which finds that the Einstein relation overcomes the effect of long temporal tails present in the pressure autocorrelation function. 30 Finally, in terms of computational effort, both methods scale linearly in time as the trapezoidal rule was used for the evaluation of the integrals in both cases, which scales with O(n) where n is the number of input data points.

C. Scaling of dynamic properties
This section analyses the viscosity values calculated from the DPD simulations and verifies that our revised Einstein formula can reproduce the rheological behaviour predicted by polymer theory. The scaling of the dynamical properties with the polymer molecular weight is used in polymer physics to determine whether the polymer solutions follow the Rouse or Zimm theoretical model. 3,51 The scaling of these quantities can be found using a linear fit and the value of the exponent is characteristic of the system's dynamic behaviour. In this work, we investigate the scaling of the specific viscosity (Z sp : where Z s = 1.1 is the solvent viscosity and Z is the viscosity of the solution) which is a unitless number that measures the polymer contribution to the solution viscosity. 47 For semidilute polymer solutions, as with the structural properties, the scaling behaviour of the viscosity depends on the values of f and f*. However, if the data are normalized by the solvent viscosity (specific viscosity) and plotted against the normalized volume fractions f/f*, the explicit effect of the polymer length  disappears and all of the data should fall onto a single master curve. 52 Fig. 5 (left) shows that indeed the specific viscosities calculated from the DPD simulations fit the same master curve. The scaling of the specific viscosity with the polymer concentration predicted by the theory is given by the following expression: and n is the Flory exponent. In athermal and ysolvent (n D 0.5), the exponent x D 2.0 and the viscosity is predicted to grow as the square of the polymer concentration.
In a good solvent where we expect hydrodynamic interactions to be dominant n D 0.6 and x D 1.25. 47 We fit the data obtained from where w E 1 indicates y-solvent while w E 0.54 good solvent conditions and D 0 refers to the diffusion coefficient at infinite dilution. In agreement with the viscosity scaling analysis, the exponents calculated from the values of D are closer to the predictions for good solvent than y-solvent conditions (Fig. 6, left).
It should be noted that in the cases with 1 o f f Ã o 2 both excluded volume and hydrodynamic interactions must be considered.
The viscosity data calculated from the melt simulations are reported in Fig. 5 (right). The friction parameter, g ij , affects the value of the viscosity and the scaling exponents. We find that the scaling of viscosity is below the expected theoretical prediction, Z B N whilst the diffusion data reported in Fig. 6 (right), follows the Rouse theory, i.e., D B N À1 the latter despite the chain lengths of some of models are quite short. 53 This difference among the viscosity and diffusion exponents has been observed previously for polymer melts modelled with similar DPD parameters as the ones employed here. In particular, Posel et al. 53 showed that the scaling of viscosity severely deviates from the scaling of the self-diffusion coefficient as the FENE spring constant increases. The authors point out that as the values of k FENE increases, the elastic contributions to the pressure tensor vanish and the viscous contributions dominate. It may therefore be expected that the scaling behaviour of the viscosities of models with high k FENE and relative short chain lengths (as those simulated here) may depart from the Rouse theory which considers only elastic contributions. Here we see the same disagreement between the D, Z exponent values, with a k FENE = 40.0 which is higher than the bond stiffness coefficients studied in the work of Posel et al. 53 However, the deviation we observe is not as severe as it is in ref. 53 Calculating the characteristic ratio, C N , (Fig. S4 in ESI †) for the melts, we also observe that, despite the higher k FENE used here, our DPD chains are more flexible than some of those modelled by Posel et al. 53,54 probably due to the different r max employed.

D. Schmidt number
The Schmidt numbers (Sc) calculated for the polymer solutions and melts using eqn (2.7) are reported in Fig. 7. It is known from the literature 55 that the Schmidt number for a DPD representation of an ideal gas (i.e., a ij = 0) can be approximated by eqn (3.6). which indicates a parabolic dependence of Sc on the friction parameter, g ij . Although we are not in an ideal gas case, this parabolic dependence seems to be followed in the present case as well, and the data points are fitted well by a function of the Bm 2 k B T where A, B are fitting coefficients (for the fittings, see Fig. S6 in ESI †). From Fig. 7, it appears evident that Sc depends on both the interaction parameters and the length of the chain, N. Irrespectively of the value of g sp , the lowest value of Sc is observed for c p = 0.3 and a sp = 25. For this conservative coefficient, the chains create ''blobs'' with a small radius of gyration (R g ) which diffuse faster within the solvent than an extended polymer chain, leading to a low Sc number. In this work, D refers to the diffusion of the chain's centre of mass rather than the diffusion of a single chain monomer. This definition is expected to increase the value of Sc compared to using the diffusion of single chain monomers. 55 Additionally, we can build a master curve for the Schmidt number plotting its normalized value (Sc/Sc 0 ) against the reduced polymer molar fraction (Fig. 8). The data can be fitted with eqn (3.7) which is obtained using the definition the   where 2 À n 3n À 1 ¼ 1:83 for good solvents and 2 À n 3n À 1 ¼ 3:0 for yconditions. As observed for the viscosity exponent coefficient, the fitting value of 2 À n 3n À 1 % 2:0, which suggests that we are closer to good solvent conditions. The term Sc 0 here is the Schmidt number calculated for infinite dilution, i.e., Finally, the ability of DPD to simulate liquid phases is evident by the fact that Sc 4 (10 2 ). which is indicative of a liquid. 12,56 The ability of DPD to simulate liquids was debated in the past 55 mainly because in DPD simulations typically Sc B O(1). Here we observe that the order of Sc is increased by a factor of 10 2 for the lowest polymer concentration, c p = 0.3 and up to 10 3 in polymer melts, c p = 1. With the correct choice of conservative interaction parameters, the modelling of high Schmidt number liquids using DPD is possible.

IV. Summary
The performance of Einstein-Helfand formalism was examined in this work in the presence of stochastic and dissipative forces and compared with the revised DPD Green Kubo (GK) formula suggested by Jung and Schmid. 18 A revised Einstein-Helfand formula was proposed for calculating the zero-shear viscosity from equilibrium DPD simulations and was applied in this work to calculate the viscosity of polymer solutions of different concentrations, c p = 0.3, 0.8, 1.0 and chain lengths N = 10, 25, 40, 50. In order to numerically verify our revised Einstein formalism a comparison with the revised GK formula suggested by Jung and Schmid 18 was carried out, and the transport properties of the systems were studied.
The viscosities were evaluated over three different conservative and friction coefficients using the time decomposition method 19 to estimate the errors. For both approaches, the viscosity was calculated as an average in the interval [2t r , 3t r ] with t r the chain relaxation time. The average viscosity value and the standard deviations from Einstein and GK were the same over 10 independent runs, however a difference of 5-8% was observed in the standard deviations for chain lengths N = 40, 50 in favour of the Einstein relation, while the average value of viscosity was the same for both methods. This difference is attributed to the error propagation in the long chains. This result indicates that the Einstein relation has higher statistical accuracy than the Green-Kubo formula at early time intervals. Additionally, we performed a scaling analysis on the polymer melts and solutions' dynamical properties and found that the scaling laws predicted by the theory are well reproduced. Our results indicate that our revised Einstein-Helfand relation can confidently be used to calculate the zero-shear viscosity for complex, slow relaxing fluids.

Data availability statement
Codes written in C++ can be found in the repository github.com/ MariaPanoukidou/zero-shear-viscosity for both methodologies applied in this work. The data that support the findings of this study are available from the corresponding author upon request.

Conflicts of interest
There are no conflicts to declare.