How does the scaling for the polymer chain in the dissipative particle dynamics hold?

We performed a series of simulations for a linear polymer chain in a solvent using dissipative particle dynamics to check the scaling relations for the end-to-end distance, radius of gyration and hydrodynamic radius in three dimensions. The polymer chains of up to 80 beads in explicit solvent of various quality are studied. To extract the scaling exponent \nu, the data are analyzed using linear fits, correction-to-scaling forms and analytical fits to the histograms of radius of gyration distribution. For certain combinations of the polymer characteristics and solvent quality, the correction-to-scaling terms are found to be essential while for the others these are negligibly small. In each particular case the final value for the exponent \nu was chosen according to the best least-squares fit. The values of \nu obtained in this way are found within the interval \nu=0.55-0.61 but are concentrated mostly around 0.59, which is very close to the best known theoretical result \nu=0.588. The existence of this interval is attributed both to the peculiarities of the method and to the moderate chain lengths being simulated. Within this shortcoming, the polymer chain in this kind of modeling is found to satisfy the scaling relations for all three radii being considered


Introduction
Variety of physical, chemical, and biological phenomena which involve polymers as well as variety of polymer characteristics one is interested in have been naturally reflected in numerous theoretical methods and models used for polymer description [ 1].In particular, it is generally recognized by now, that the scaling properties of a long flexible polymer chain immersed in a good solvent are perfectly described by the model of self-avoiding walks (SAW) [ 2].A textbook example is given by the mean square end-to-end distance R 1N of a SAW of N steps and that of a linear polymer chain, either of which scales in an asymptotic limit N → ∞ as: with a universal exponent ν.The value of the exponent depends on the space dimension only and is the same for all SAW on different three dimensional (3d) lattices.Relation (1) is violated for the finite N and an approach to the asymptotics is governed by the corrections to scaling.With an account of the first correction equation (1) reads: Here A, B are non-universal amplitudes and ∆ is an universal correction-to-scaling exponent.Theoretical estimates for the exponents follow from the field theoretical renormalization group calculations [ 3]: ν(3d) = 0.5882 ± 0.0011, ∆(3d) = 0.478 ± 0.010.
c J.M.Ilnytskyi, Yu.Holovatch Whereas the SAW on the discrete lattice is a perfect model to study scaling of long flexible polymer chains in a good solvent [ 4], for obvious reasons it fails to describe most of the other properties of the chain [ 1,5].A more realistic model would be defined in a continuous space and involves both monomer-monomer and monomer-solvent interactions.Brownian dynamics is one of the candidates, as it operates on the mesoscale and incorporates the random force acting on each monomer which mimics the effect of a solvent [ 6].However, the method lacks correct hydrodynamics limit.Similar but more rigorous technique, i.e., the dissipative particles dynamics (DPD), was introduced by Hoogerbrugge and Koelman [ 7] and later refined by Español and Warren to satisfy the detailed balance [ 8].The method has correct hydrodynamics limit [ 9] and quickly gained much attention as a powerful technique for mesoscale simulations of various kinds of macromolecules [ 6].
Therefore, an important issue of analysis is to check how the scaling laws (e.g. of equation ( 1) form) hold in the case of the DPD-based simulations.This issue was partly considered in [ 10,11,12,13,14].However, some issues remain unclear (see the next section 2 for details).We see a good reason to discuss in detail the scaling laws for the polymer chain in a solvent in the DPD method.As far as the method was successfully used for the simulations of branched, amphiphilic and other complex molecules, the validation of correct scaling laws on such a well-studied system as linear chain polymer would be a strong supporting factor for extending the subject to the scaling laws in more complex systems (e.g.star-like, branched, dendritic systems, etc.).
The set-up of our paper is as follows: in the next section 2 we give a brief review of former simulation results, section 3 contains the details of our simulations.Our main results are presented in section 4, section 5 gives conclusions and outlook.

Previous simulations
A number of various simulation techniques have been used to numerically verify the scaling laws of the polymer in solution, namely the molecular dynamics (MD), Monte Carlo (MC), DPD and others.Here we shall provide a rather brief account of references relevant to the present study.
By using the MC method enriched by special techniques, one can move very efficiently through the conformational space (see, e.g.[ 6]).As the result, the method is quite suitable for the simulations of random walks (RW) and SAW [ 15,16].Li et al. [ 16] used a pivot algorithm and other techniques applied to the chain of N = 100 ÷ 1000 monomers.The value for the scaling exponent obtained is ν = 0.588 and it agrees very well with the theoretical estimate (3).To study dynamical properties, however, one should turn to true dynamical methods.Pierleone and Ryckaert [ 17] performed a comprehensive analysis of the scaling laws for both static (radius of gyration, R g and end-to-end distance R 1N ) and dynamical (diffusion coefficient, D) properties using MD simulations.One of the important outcomes is that the static properties are unaffected by the finiteness of the box size L providing that L/R g > 3. The scaling exponent values ν = 0.59 and 0.568 (for the R 1N and R g , respectively) were obtained.On the contrary, dynamical properties were found to be strongly dependent on L. This fact was explained by hydrodynamic interactions of the chain with its own images due to the periodic boundary conditions (PBC).If, following Dunweg and Kremer [ 18] one uses the Kirkwood formula, then a correct dynamical scaling is also recovered [ 17].
A number of simulations of linear polymers using DPD technique are available.Schlijper et al. [ 10] were the first to study static and dynamic scaling laws for the DPD chain in athermal solvent.Later their results for this case were confirmed and extended to include the solvents of variable quality by Kong et al. [ 11].These authors studied R g (among the other properties) and have found for the athermal solvent ν = 0.52 which is closer to the polymer at the θ-condition.With the increase of solvent quality, the value ν ≈ 0.6 close to a SAW exponent was recovered.This problem was further addressed by Spenley [ 12] who studied the R 1N in polymer melt and solution.In the latter case, relevant to our study, the exponent ν = 0.58 was obtained.The previous result ν = 0.52 by Kong et al. [ 11] was commented as being influenced by the type of spring potential used there.In our opinion, this discrepancy could also be due to insufficient statistics.
It has also been debated whether the soft repulsion inherent to the DPD method is sufficient to model the self-avoidance of the chains.In particular, Symeonidis et al. [ 13] used additional intramolecular interactions (e.g.Lennard-Jones-like repulsion term) to increase the chain selfavoidance and experimented with different forms of the bond springs (i.e.FENE, Hookean and WLC).One of the conclusions made by the authors is, that the presence of FENE potential alone ensures a correct value for the exponent ν ≈ 0.6 and no extra interactions for the improvement of self-avoidance are required.Jiang et al. [ 14] have carried out an extensive study of the hydrodynamical properties obtained via the DPD method.In particular, scaling laws for both R g and R 1N for the dilute athermal solution are relevant to our study.The authors found the values for the scaling exponent ν weakly dependent on the box size with the averages of ν = 0.58 for the radius of gyration and ν = 0.595 for the end-to-end distance [ 14] (the box sizes of L = 15÷30 DPD length units were used and a range of chain lengths was N = 10 ÷ 100 DPD beads).The other conclusion given in [ 14] is the one already mentioned above, claiming that the original soft interactions are sufficient for self-avoidance of the chains in the DPD method.Similarly to the MD study of Pierleone and Ryckaert [ 17], the diffusion coefficient was found to be highly sensitive to the simulational box size.However, after the appropriate corrections to the system size were introduced, the scaling law for the diffusion coefficient D ∼ N −ν was obtained.Since, according to Kirkwood theory, the diffusion coefficient D is related to the hydrodynamic radius R h [ 18], the latter is also shown to obey the correct scaling law.
Despite the availability of simulational studies that concentrate on static and dynamical scaling of the polymer in DPD method, there are still certain points that need to be clarified, namely: (i) do all the relevant radii, R 1N , R g and R h scale with the same scaling exponent ν?
(ii) are the corrections to scaling relevant and can these improve the results or make these more consistent?
(iii) in case of a good solvent, will the changes in the polymer-solvent interaction potential lead to any detectable drift of the scaling exponent ν?
These three points will be addressed in the current study.

Simulational details
In our study we closely follow the DPD method as discussed by Groot and Warren [ 19].Polymer molecule (see, figure 1) is modeled as a chain of soft beads connected via harmonic springs.The bonding force acting on the i-th bead from its bond neighbour j is: where Here and thereafter, the length, mass, time and energy (k B T ) units are chosen to be equal to unity.The solvent is described in an explicit way, in a form of isolated beads of the same origin as the polymer ones (see, figure 1).Each i-th bead (polymer or solvent one) is subject to a pairwise non-bonded force from its j-th counterpart.The force acting on the i-th bead due to the interaction with the j-th bead consists of three contributions: where the conservative F C ij , dissipative F D ij and random F R ij contributions are of the following form: Here v ij = v i − v j , v i and v j being the velocities of the beads, θ ij is Gaussian random variable: According to Español and Warren [ 8], the dissipative and random force amplitudes are interrelated, w D (r ij ) = (w R (r ij )) 2 and σ 2 = 2γ, to satisfy the detailed balance.The frequently used analytical form is: We choose the value γ = 6.75.Parameter a in the conservative force F C ij defines maximum repulsion between two beads which occurs at complete overlap r ij = 0. Three different a parameters can be distinguished, a pp for the polymer-polymer, a ss for the solvent-solvent and a ps for the polymersolvent interactions and all these can be varied.As was shown in [ 19], these parameters can be related to the Flory-Huggins parameter χ.To integrate the equations of motion, we used the modified velocity-Verlet algorithm by Groot and Warren [ 19].Let us note that the structural properties of the chain in a solvent are defined solely by the conservative forces (6), while the dissipative and random forces (7) govern the phase space trajectory.

(9)
The terms "good" and "very good" are indicative rather than exact, since only repulsive forces are taken into account in DPD.Hence, "good solvent" means less repulsion for the polymer-solvent interaction than that for the polymer-polymer one.The last two cases in equations ( 9) correspond to the values ξ = −0.2 and ξ = −0.6,respectively, in notation of Kong et al. [ 11].
All the simulations were performed in N V T ensemble.At first, the initial chain conformation was generated as a RW.Then, it was immersed into a relatively large box filled by solvent and short preliminary runs were performed to estimate the approximate R g for the polymer.Thereafter, the final run was performed using the cubic box of the linear size L ≈ 5R g .This size is considered to be sufficiently large to eliminate the effects of PBC [ 14,17].The number of solvent beads N s was chosen to keep the total reduced density of the solution equal to ρ * = (N + N s )/V = 3.The time step was chosen to be ∆t = 0.04 in reduced units.We performed from 3 • 10 6 (for N = 5) up to 15 • 10 6 (for N = 80) DPD steps for the final run.We did not concentrate on the relaxation times in detail, but found the behaviour similar to Jiang et al. [ 14], τ ≈ 0.13 N 1.81 who used similar DPD parameters for athermal solvent.This brings us to the following estimates for the simulation runs: about 5 • 10 4 τ for N = 5 and 1.7 • 10 3 τ for N = 80 with the estimates for intermediate values of N falling in between.The chain conformations were saved after performing every 1000 DPD steps and last 3/4 of these were used for the subsequent analysis.As far as all the simulations were performed using PBC, the unwrapping of the chain at each time step was performed based on the connectivity information.

Results
We concentrated our analysis on the three following metric characteristics of the linear polymer.The instantaneous squared end-to-end distance at a given time t is defined as: where r 1 and r N are the radius-vectors for the first and last bead of the chain, respectively.The squared radius of gyration at time t is defined as: where R COM is the radius-vector of the polymer center of mass.The inverse hydrodynamic radius is defined as follows: where the sum runs over all different bead pairs i, j.All three properties are time averaged (denoted hereafter as • • • ) and then the first two can be square rooted while the third one is inverted to provide the estimates for R 1N , R g and R h , respectively.One should consider R g and R h as more useful characteristics as compared to R 1N .First of all, R 1N is meaningful for the linear polymer or for a certain subchain of the branched molecule (e.g.backbone, side chain, star-polymer arm, etc.), while both R g and R h are universally defined.Apart from that, both these properties are directly measurable via small-angle scattering experiments and via small-angle dynamic light scattering experiments, respectively (for more discussion on this, see [ 20]).

Linear fit and first correction-to-scaling fit
As was already mentioned above, a number of universal scaling laws have been observed for the polymer properties [ 2] in the limit of infinite polymer length, e.g.equations (1).To reproduce these properties, one needs to consider rather long chains in the simulations [ 16].Due to the polymer coil self-similarity, one would expect that such a coarse-grained approach as DPD should be capable of reproducing correct scaling laws for much shorter chain lengths.For the moderate polymer lengths, usually the scaling is written in terms of a number of bonds, N − 1, instead of the number of monomers N .Also, as is known, the correction-to-scaling terms are important for moderate system sizes, as is indicated by the studies of the phase transitions in spin lattice models [ 21,22,23].Taking into account the first correction-to-scaling term, equation ( 2), one has for the We shall work with square rooted values and keep the linear term only in the correction-to-scaling expansion above.Taking a logarithm of both sides of equations ( 14), we obtain ln where A ′ and B ′ are self-explanatory.The same formulae can also be applied for R g = R 2 g and R h = 1/ 1/R h (in the latter case one has B ′ ≈ −B).As a result, the same fit (15) can be used for all three quantities under interest, R 1N , R g and R h .In the case of linear fit the correction-toscaling amplitude B ′ is assumed to be zero and one recovers the linear fit form.To perform the fits we used the least-squares routine mrqmin from the Numerical recipes book [ 24].The accuracy of the fit was monitored by the cumulative squared deviation per one data point, where (N [k] ,R 1N [k] ) is the k-th out of n data data points and F (N ) is the fitting function.It is generally known that very good statistics for each data point are required to ensure robust results for the scaling exponent ν.Otherwise the latter depends essentially on the selection of data points being used for the fit.This was also found to be the case in the course of the present study.Hence, relatively lengthly simulations are performed to improve statistics.To check the consistency of the results, we performed a series of fits using a gradually increasing range of intervals in the chain length N , e.g.N = 5 ÷ 28, N = 5 ÷ 40, • • • N = 5 ÷ 80, which are presented below.
As was anticipated, both cases of "athermal 1" and "athermal 2" solvents were found to show the same averaged properties within the accuracy of the simulations.Therefore, the data for both cases have been averaged at each N and marked as "athermal" further in the text.At first glance, the data points for both R 1N and R g are well fitted by the line in log-log scale (see, figure 2), but, in fact, the best fits were achieved using the first correction-to-scaling form (15).We should remark here that we were unable to achieve stable results for both R 1N and R g using the four-parameter fits, in which all A ′ , ν, B ′ and ∆ in (15) are fitted simultaneously.Instead, we were forced to fix the correction-to-scaling exponent ∆ at certain value and to perform the fit over three remaining parameters.We used two choices for ∆, namely the best theoretical value ∆ = 0.478 [ 3] and ∆ = 1.The results presented in table 1 indicate that the correction-to-scaling term for R 1N is essential in the case of athermal solvent only.In the table, we also show the linear fit performed over the largest chain lengths, N = 20 ÷ 80, for the sake of comparison.One can see that the latter leads to almost the same result ν = 0.586 as the correction-to-scaling fit ν = 0.591, which is indicative of the crossover quickly decreasing at N > 20.In the cases of good and very good solvents very small amplitudes B ′ are found and, as the result, linear fits are almost equally good (see, table 1).For R g the situation is exactly reversed.The correction to scaling is important for the cases of good and very good solvents (see, table 2).We cannot suggest some definite explanation for this effect, except the conclusion that preference should be given to the general form (15), rather than to the linear fit, otherwise the exponent ν could be essentially under-or overestimated (see, tables 1, 2).The four-parameter fit, mentioned above, does not work for either R 1N or R g but is supposed to perform better for R h , where essentially larger correction to scaling is to be expected [ 20].Indeed, this is found to be the case in our simulations.First of all, this can be noticed from the arrangement of the data points, see figure 3. The results based on numerical fits are presented in table 3, where both the cases of fixed ∆ = 0.478 and of the four-parameter fits are shown.The latter give essentially better results in terms of fitting accuracy and provide independent estimates for both exponents, ν ≈ 0.56 ÷ 0.58 and ∆ ≈ 0.65 ÷ 0.70.As was demonstrated by Dunweg et al. [ 20] there are two comparable correction-to-scaling terms for the R h , one is proportional to N −∆ and another, the so-called "analytic" correction, is proportional to N − (1−ν) .Therefore, if one uses the form (15) with a single correction term, then a rather effective exponent ∆ will be obtained.This is how we should interpret the value ∆ obtained in our study.

The fits based on the distribution of R g values
The averaging of the metric properties of interest, e.g.R g = R 2 g , throughout this study was performed as the arithmetic mean over discrete set of time points.One can also be interested in the distribution p(R g ) of the R g values themselves and, for instance, in examination of the scaling law for the most probable values with the chain length N .The following analytical form has been postulated by Lhuillier in 3d [ 25] p where α = 1/(3ν − 1), δ = 1/(1 − ν) and A 1 , A 2 are constants.Subsequently, this has been verified in a number of lattice MC simulations [ 26,27].One can rewrite this distribution in terms of R 2 g , using the expressions for α and δ and incorporating the N -dependent factors into the constants, Table 2. Linear fit and first correction-to-scaling fits to the form (15) for Rg with fixed exponent ∆ = 0.478 and ∆ = 1 when using various sets of data points.The best fit for the interval N = 5 ÷ 80 is underlined.
fitting linear fit fit with ∆ = 0.478  and the result reads: This distribution has a maximum at which can be square rooted to obtain an estimate for the most probable value of the radius of gyration R g max = [R 2 g ] max .We prefer to work with the distribution of R 2 g (18) and then to square root the maximum to match the averaging procedure for simple averaging, R g = R 2 g .From the point of view of numerical fitting, the form (18) has four parameters: one exponent ν and three coefficients, a 1 , a 2 and C. Therefore, one can have a first estimate for the exponent ν right from the fit of the data to the form (18).The second estimate for ν can be obtained from the anticipated scaling of R g max according to the form (15).The distributions p(R 2 g ) were built in a form of histograms and were found to map very well on the Lhuillier form (18), see samples in figure 4 obtained for the case of athermal solvent.These fits, as was mentioned above, provide both the estimate for R g max and an independent estimate for the exponent ν.The latter was found to be scattered in the interval ν = 0.53 ÷ 0.58 depending on N and solvent type.We do not expect high accuracy from such a fitting, as far as dependence on ν in equations ( 18) is rather complex.Nevertheless, the values of ν found are very reasonable and indicate a self-consistency of this numerical approach.The next step is to use the maxima positions R g max of the distributions (18), found via equations (19), and to find a fit to the scaling law (15) similarly to the analysis of R g data in the previous subsection.The data points alongside with the best fits are shown in figure 5 and the results for the exponent ν, correction-to-scaling amplitude B ′ and fitting error χ 2 are presented in table 4. We would like to stress that, contrary to the results given in table 2, no dependence on solvent quality is observed in table 4 and the same value ν = 0.582 − 0.585 is found essentially in all three cases.This might indicate that the effect of small drift of ν observed in a very good solvent in table 2 could be an artefact of the numerical method used.The cumulative outcome of our analysis for the scaling exponent ν which governs the scaling law (2) is presented in figure 6.Here we show a histogram of the values of ν obtained for R 1N , R g and R h of a polymer chain in different solvents considered in the current study, namely, in athermal, good and very good solvents.To make this histogram, the best fits (underlined results in tables 1-4) have been used.One can see from the histogram that the value of ν is scattered in an interval ν = 0.55÷0.61but in fact is found predominantly in much narrower interval ν = 0.58÷0.60centered around the best theoretical estimate ν = 0.5882 ± 0.0011 [ 3].This allows us to conclude that within the accuracy of our analysis, the scaling law (2) holds for all polymer characteristics and in all solvents considered in this study.Our conservative estimate for the scaling exponent therefore is ν ≈ 0.59 ± 0.01. (20)

Conclusions and outlook
In this study we performed DPD simulations of the polymer chain in a solvent of various quality with the aim to reexamine how well the scaling laws hold for various polymer characteristics.Chains of up to 80 beads were considered and simulation runs from 1600 up to 50000 relaxation times were performed.Three metric properties were considered, end-to-end distance, radius of gyration and hydrodynamic radius.For the analysis we used the linear fits, first correction-to-scaling fits and fitting the distribution of the radius of gyration.As an outcome, within the accuracy of the simulations and of the data analysis technique, the following conclusions can be drawn: (i) All three metric properties obey the scaling law with very close values for ν found within the interval 0.55 ÷ 0.61 (depending on the method of analysis).Most values are concentrated around the average ν = 0.59 being very close to the theoretical estimate ν = 0.5882 ± 0.0011 [ 3], see, figure 6.
(ii) Corrections to scaling are found to be important in all combinations except a few property/solvent ones and thus the conclusion (i) is valid only when the correction-to-scaling terms are taken into account.
(iii) No or very small (up to 4%) drift of the scaling exponent ν is observed when changing the solvent quality (remaining, however, in a good solvent regime), but this effect depends on a numerical technique being used.
The method of analysis can be extended to the study of the scaling laws in more complex molecular architectures, e.g.star-like polymers, branched and hyperbranched molecules (including amphiphiles).

Figure 1 .
Figure 1.The polymer molecule in a solvent, as modeled in the DPD method.Both monomers constituting a polymer chain and solvent molecules are considered as interacting soft beads (grey and hollow discs, correspondingly).

Figure 2 .
Figure 2. Data points for R1N (on the left) and Rg (on the right) and their best fits using the form (15) obtained from the DPD simulations on a set of chain lengths N = 5 ÷ 80 at different solvent conditions, see notations on both plots.

Figure 3 .
Figure 3. Data points for R h and their best fits using the form (15) obtained from the DPD simulations on a set of chain lengths N = 5 ÷ 80 at different solvent conditions, see notations on the plot.

Figure 4 .
Figure 4. Mapping of the distribution of radius of gyration values p(R 2 g ) obtained from the DPD simulations in athermal solvent, onto Lhuillier form (18), a few indicated chain lengths are shown only for the sake of clarity.

Figure 5 .
Figure 5. Data points for Rg max and their best fits using the form (15) obtained from the DPD simulations on a set of chain lengths N = 5 ÷ 80 at different solvent conditions, see notations on the plot.

Figure 6 .
Figure 6.A histogram of values of ν obtained for different characteristic radii of a polymer chain in different solvents(9).The data are collected according to best fits results (underlined in tables 1-4).The value of ν is found predominantly in the interval ν = 0.55 ÷ 0.61 and concentrate around the best theoretical estimate ν = 0.5882 ± 0.0011[ 3].

Table 1 .
Linear fit and first correction-to-scaling fits to the form (15) for R1N with fixed exponent ∆ = 0.478 and ∆ = 1 when using various sets of data points.The best fit for the interval N = 5 ÷ 80 is underlined.

Table 3 .
Linear fit and first correction-to-scaling fits to the form (15) for R h with fixed exponent ∆ = 0.478 and from the four-parameter fits when using various sets of data points.The best fit for the interval N = 5 ÷ 80 is underlined.

Table 4 .
(15)ar fit and first correction-to-scaling fits to the form(15)for R max