Roland Leonel

Nkepsu Mbitou

Alain Dequidt

Florent Goujon

Benoit Latour

Julien Devémy

Nicolas Martzel

Patrice Hauret

Patrice Malfreyt

Development of anisotropic force fields for homopolymer melts at the mesoscale

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

I. INTRODUCTION

Coarse-grain (CG) models have been used for a long time for the molecular simulation of polymers, in order to reach longer time scales and bigger systems [1][START_REF] Schweizer | Rajasekaran[END_REF][3][START_REF] Kremer | Computer Simulation in Chemical Physics[END_REF][START_REF] Akkermans | [END_REF][6][7][8][9][10][11][12] . The first CG models were quite rudimentary and used to focus on generic properties of polymers, like scaling laws 10,11,13,[START_REF] Rapaport | The Art of Molecular Dynamics Simulation[END_REF] . Nowadays, the standards have raised and while researches are still active on the generic properties of polymers, there is a growing demand for simulations of specific polymers with quantitative predictions of thermo-mechanical properties. The challenge is to develop force fields at the coarse-graining level, that are at the same time fast to simulate and as accurate as possible in reproducing some experimental properties. [START_REF] Gkeka | [END_REF][16][17] The most realistic CG force fields (FF) distinguish between grain types and use tabulated effective interaction potentials, because they form an expressive family of functions. They also consider random forces and a viscous friction to account for the contribution of the omitted fast degrees of freedom at small scale [START_REF] Grabert | Projection Operator Techniques in Nonequilibrium Statistical Mechanics[END_REF]19 . The shape of the grains is also a key factor. At order 0, grains are point particles located at the center of mass of the underlying atoms, and their isotropic interaction potential fixes their size. Most of the CG models use this level of approximation. At order 1, grains are ellipsoids whose shape corresponds to the gyration tensor of the underlying atoms. The theoretical interest for such models has existed for a long time, but until now in general, the gain in realism has not justified the extra complexity Historically, anisotropic models have been widely used to simulate liquid crystals, where they are essential [20][START_REF] Faller | Reviews in Computational Chemistry[END_REF][START_REF] Wilson | [END_REF] , but a few applications to polymers exist [23][24][25][26] . For even more realism, it would be possible to develop more complicated but more accurate grain shapes 27 at the expense of more computation time. One could argue, that by splitting a) Electronic mail: alain.dequidt@uca.fr b) Electronic mail: florent.goujon@uca.fr the grain into smaller grains one obtains a more realistic shape, which is generally true, but the computation time would increase and the benefit of coarse graining would disappear. Another option is to use machine learning to build many-body anisotropic potentials 28 .

In this paper, we choose to use a very simple anisotropic model, where each grain is an ellipsoid of revolution (spheroid) of fixed shape, whose main axis is oriented parallel to the polymer backbone 29 . This model has the advantage of simplicity and efficiency. Of course, this is still a crude approximation and more realistic ellipsoidal models are under development, including shape fluctuations and general orientations 30 .

Bigger grains are interesting because they reduce more the number of degrees of freedom. At the same time, the effective interaction potentials become softer because big grains are allowed to overlap without implying that the underlying atoms overlap. This has major consequences: 1) The CG models are usually much more compressible than the fine grain models. 2) More dramatically, the polymer chains may cross each others, preventing efficient entanglements to form, whose impact on the chain dynamics are fundamental.

By using more realistic grain shapes, we expect that: 1) The resulting CG model accounts for the chemical specificity of the polymer 2) Even with strong potential, the grains can get close by taking suitable orientations. Compressibility could therefore decrease by using ellipsoidal grains. 3) The elongated grain shape would "fill the gap" between consecutive grains along the polymer backbone, hindering the nonphysical chain crossing. This effect was demonstrated using a generic model in ref. 29.

The purpose of this work is to investigate to what extent these expectations are true. What is the impact of grain anisotropy on the model realism at different coarsegraining levels?

In the following, we studied the effect of non-spherical shapes of coarse-grained (CG) particles, on static and dynamic properties of a 1,4-cis-polybutadiene (cPB) system melt. Several CG models of cPB can be found in the recent literature 25,[31][32][33][34] . Most often, the CG potentials are obtained using Iterative Boltzmann Inversion (IBI) 35 , which aims at reproducing the local structure of the underlying atomistic systems. By using this method, the local structure is therefore almost perfectly reproduced. With IBI, obtaining the correct density requires to tune a "pressure correction" and is not straightforward. On the other hand, compressibility is not tuned explicitly. It is widely recognized 33,34,36 that CG models are too "soft", above all at high coarse-graining level, which leads to too a high compressibility. A possible remedy would be to use local-density-dependent potential 37,38 , but these are hard to parameterize, especially in mixtures. Usually, the compressibility issue is simply disregarded or simulations are run in the N V T ensemble. Such an approach is not satisfactory for studying mechanical properties and may result in bad local density, especially with mixtures. Dynamical properties are generally obtained by artificially scaling the time of CG simulations 31,34,39 . However, as pointed in Ref 34, there is no guarantee that "a single scaling would reproduce all dynamical properties". Ref 32 is more original. Here, the dynamical properties are obtained using the method of Hijòn et al. 19 , based on force fluctuations. The method is tested on small, unentangled chains only (12 monomers), for which bond crossing is not a dramatic issue. It requires to run constrained simulations at the atomistic level. This reference (as well as Ref 40) points out that the grain anisotropy should probably be considered for more consistency and has motivated the work of Ref 25. The latter reference introduces an anisotropic CG model of cPB. This model is derived using a method similar to Ref 32 and requires to run a constrained dynamics. It makes use of torques and angular velocities. The chains are also small and not entangled, so that potential bond crossings are not an important issue. The model is validated based on structure only and does not consider pressure or density at the CG level. Finally, several methods have been proposed in the literature to deal with the unphysical bond crossings that are observed at high CG level and for long polymer chains 29,34,[41][42][43][44][45] . These methods are most often applied to generic systems, but sometimes a quantitative agreement with a specific system is researched, using multiscale approaches 34 .

In the present study, we use the backbone oriented anisotropic model of Ref 29 to assess its potential in producing realistic CG models of possibly entangled cPB. The choice of this model is motivated by its simplicity and numerical efficiency in line with the mesoscopic simulation method used here, and because we have at our disposal a bottom-up method for parameterizing it optimally, namely Statistical Trajectory Matching 46 .

In the first section, we introduce the model and explain how the force field parameters are obtained. In the second section, we present the results and discuss the effectiveness of the model.

II. METHODS

A. Coarse grain model

The anisotropic grain model used in this work is the so-called backbone-oriented anisotropic model (BOA) described in ref. 29. Here we briefly recall its principle.

Since cPB is a homopolymer, only one kind of grain is defined. Each grain is represented by an ellipsoid of revolution of fixed shape. The aspect ratio of the grain is denoted by h. It is natural to ask if the value of h can be obtained by directly analyzing the gyration tensor of the grain in MD simulations. This is not so easy, because the actual grain shape is not constant in time, especially for multi-monomer grains, and does not actually have a revolution symmetry 40 . In this study, we have tested various values of h between 1 (isotropic grains) and 4.

The axis of grain i is ⃗ u i ∥ (⃗ r i+1 -⃗ r i-1 ), where ⃗ r i is the position of grain i and i -1, i + 1 are the neighbor grains along the polymer chain.

Bond interactions are harmonic and isotropic. Angular interactions are also harmonic (see Supporting Information). Non-bonded interactions are pairwise and anisotropic: The pairwise potential is a tabulated function U (R) of the effective distance R between the grains. This orientation-dependent effective distance is defined for grains i and j by

⃗ R ij = h 1 6 I 3 + (h -1) ⃗ u i ⊗ ⃗ u i + ⃗ u j ⊗ ⃗ u j 2 -1 2 • ⃗ r ij (1)
where I 3 is the identity matrix. When h = 1, ⃗ R ij ≡ ⃗ r ij is the true distance between the grain centers. Otherwise, when the two grains are parallel side by side, their effective distance is greater (R > r), while when they are parallel in single file, their effective distance is lower (R < r). For the computation of forces, the grain orientations are assumed to be fixed so that no torque is considered. During such computations, the force is interpolated using cubic splines through 10 points at fixed distances within the selected cutoff radius (see Supporting Information).

Thermalization and the effect of fast degrees of freedom on the dynamics are obtained thanks to a Dissipative Particle Dynamics (DPD) thermostat. The parameters of the tabulated non-bonded force, of the bonded interactions and the friction are obtained using the bottom-up method called Statistical Trajectory Matching (STM) 40,[46][47][48] . This is the first time this method is applied to the BOA model.

B. Atomistic simulations

The bottom-up method requires high resolution trajectories obtained at the atomistic level. Here we give the simulation details of these atomistic simulations. The generated cPB system is a cubic simulation box, and it consists of one periodic amorphous cell. Specifically for the atomistic trajectories necessary for the STM process, such system contains 90 all-atom (AA) cPB chains with 60 monomers of 1,4-cis-butadiene per chain. Each chain is successively generated in the initial configuration using a homemade Python code, by imposing the bond and the angles parameters. The resulting amorphous system is finally accepted if its density is equal to the imposed value, and the averaged end-to-end distance is close to the experimental one 49 . Molecular Dynamics (MD) simulations are performed through the LAMMPS software 50 , and the generalized Amber force field 51 (GAFF) is exploited to model the polymer. The parameters of this force field (FF) are given in the Supporting Information. The simulation protocol is elaborated by first making an equilibration phase of 10 ns. Secondly, a short production phase of 1 ns is performed to record configurations of the atomistic trajectory every 25 time steps, by using a time step of 2 fs. During these phases, the resolution of the fundamental equations of motions is carried out by using the standard velocity Verlet algorithm 52 . The oscillations of all the C-H bonds are constraint in the polymer system, by using the SHAKE algorithm 53 . MD simulations are performed in the constant N pT ensemble, where N , p, and T are the number of atoms, the pressure and the temperature, respectively. The atmospheric conditions (p = 0.1 MPa, and T = 300 K) are then maintained thanks to the Nosé-Hoover barostat and thermostat algorithm 54,55 , where the barostat and thermostat relaxation times are set at 2000 and 200 fs, respectively. The cutoff radius of the Lennard-Jones (LJ) and Coulombic (at long-range) interactions is set at 12 Å, and the (1-2), (1-3), and (1-4) weighting factors are set at 0.0, 0.0 and 0.5, respectively. For the LJ interactions, the crossing terms of different atoms i and j are computed using the Lorentz-Berthelot mixing rules as follow

ϵ ij = √ ϵ i ϵ j and σ ij = 1 2 (σ i + σ j ), (2)
where ϵ ij and σ ij represent the energy parameter and the diameter of i and j atom types, respectively. At a long-range, the Coulombic interactions are handled by using the particle-particle-particle-mesh (PPPM) 3D method 56 . Finally, the periodic boundary conditions are applied in all the three directions of the system. Configurations of the obtained atomistic trajectory are recorded every 50 fs during a short MD production phase, and such time is chosen as the time step of the next CG simulations. In order to get different high resolution trajectories, each polymer chain of the atomistic configurations is mapped into 5 different coarse-graining levels λ (where λ varies from 1 to 5, in steps of 1). The value of λ represents the number of repeat units, which are mapped into one spherical CG bead located at the center of mass of all the constitutive atoms. The figure 1 illustrates the mapping of one cPB chain in three differ- ent coarse-graining levels.

C. Parameterization of the mesoscopic interactions

The parameters of the CG model are obtained using the STM method 40,[46][47][48] . The purpose of the article is not to promote this method with respect to al- ternative older methods (iterative Boltzmann inversion, force matching,...). We note however that applying this method to anisotropic interactions is straightforward.

Here we briefly recall how STM works.

Any set of parameters gives a CG model, using which CG trajectories can be generated. Because of the random forces of the DPD model, a whole set of trajectories could be generated starting from any initial configuration. Among all these trajectories are the reference trajectories obtained from high resolution simulations. By assuming Gaussian random forces, the probability for a given model to generate exactly the reference trajectories can be computed analytically. The STM method gives the model parameters, using which this probability is the highest.

In practice, instead of maximizing the probability P, we minimize

L = ⟨-ln P t ⟩ t (3)
where the notation ⟨⟩ t means the average over time steps in the reference trajectory and P t is the probability to generate the next configuration exactly as in the reference trajectory.

The range of the non-bonded interaction is a non-linear parameter which is not optimized by the STM method. We have tested several values of the cutoff distance r c between 20 and 40 Å and finally retained r c = 20 Å. On Figure 2.c, we show that this choice is effectively the best for λ = 3. We also note that the most realistic value of h in the sense of a better likelihood is typically between 1 and 2.

D. Coarse grain simulations

The STM method is applied to find the optimal parameters of the CG FF. Then, 65 CG models are obtained separately by varying 13 values of h (ranging from 1 to 4 in steps of 0.25), and 5 values of λ (ranging from 1 to 5 in steps of 1). The procedure of the previous section is used to generate one AA cPB system, which contains 20 chains with 900 monomers per chain. In accordance to the value of λ, such cubic and amorphous system undergoes a mapping process to get initial configurations at the mesoscale. DPD simulations 57 are carried out at this scale, using the LAMMPS software. The simulation protocol begins by an equilibration stage of 10 ns, it is followed by a short production stage of 1 ns (performed to record configurations of CG trajectories every 1 time step), and it ends by a long production stage of 1000 ns (in this stage, configurations are recorded every 2000 time steps). A time step of 50 fs is used during these simulations, and they are performed in the constant N pT ensemble. The atmospheric conditions (p = 0.1 MPa, and T = 300 K) are then maintained by using both the Berendsen barostat algorithm 58 (where the bulk modulus and the damping parameter are set at 100 MPa and 50 ps, respectively), and the DPD thermostat. For conservative interactions, the cutoff radius is set at 20 Å, whatever the coarse-graining level λ. Moreover, the weighting factors for pairwise conservative energies and forces between bond and angle patterns are set at 0.0 and 0.5, respectively . For non-conservative interactions, the cutoff radius is also set at 20 Å, except for λ = 1 which is defined at 10 Å. Consequently, DPD forces act on each CG bead of any system, and they are integrated by using a standard velocity Verlet algorithm.

III. RESULTS AND DISCUSSIONS

A. Coarse-grained potentials

The CG effective potential was optimized for each coarse-graining level λ and anisotropy factor h. The re- sult is shown in Figure 3. Note that in these figures, the energy is plotted as a function of the effective distance R, not r.

It can be seen that the usual potential shape is obtained in all cases, with repulsive interactions at short distance and a single smooth attractive well at longer distance. The global trend is that the more h increases, the more the potential depth decreases and the more the potential well shifts to smaller R. The potential depth decreases by more than a factor 2, which is really significant. It means that attraction need not be as high at long distance as in isotropic models in order to keep the system density at constant pressure.

The excessive softness of the coarse grain model is more related to the repulsive interaction at short distance. The comparison is difficult because R depends on the grain orientations. Note that the softness of the interaction is not trivially related to position of the well. In fact it can be observed in Figure 4 that the nearest neighbors are around r = 6 Å to 7 Å, namely in the repulsive range, not close to the minimum of the potential. We believe, that this is a many-body effect: the weak attraction at longrange pushes the nearest neighbors close to each other, at a distance where they (softly) repel.

The hardness of the interaction may be defined as the stiffness, ie the second derivative of the potential, at the average distance of the closest neighbors (first peak of the radial distribution function). This hardness is not directly related to the width of the well. This definition may be effective for isotropic potentials, but here not only R depends on the grain orientations, but so do the derivatives of U with respect to r. As a consequence we are not going to establish here the effect of anisotropy on the softness or hardness of the potential. Instead we delay this discussion to the section about compressibility.

B. Structures

The impact of the anisotropy factor h on the structure at different λ is shown on Figure 4 by means of the (isotropic) radial distribution function (RDF). As expected, the distance of the nearest neighbors increases with λ due to the bigger size of the grains. The structure also becomes less sharp and the peak of nearest neighbors is less high. This is also expected and related to the increasing "softness" of the interactions. It can be observed, that the quality of the structure reproduction does depend on h. At λ = 1, the RDF is very well matched using h ≈ 1.5. At greater λ, the RDF is still well reproduced, but the agreement is less good. For a quantitative measure of the structure matching, we computed the discrepancy defined by err =

20 0 |g(r) -g ref (r)| dr (4)
This is shown in Figure 5. Following this criterion the optimum is h = 1.75 (resp. h = 3.25, and h = 4) at λ = 1 (resp. λ = 2 and λ = 3).

We emphasize that the STM method is not based on matching radial distribution functions. The agreement is just used to estimate the quality of the CG force field in reproducing the structure a posteriori.

Although the bonded interactions remains isotropic in this model, h may also impact the distribution of bonds, angles. . . The distributions of bond distances and angles between consecutive segments is given in the Supporting Information. From the corresponding errors, it is found that the higher h, the better the bond distribution, but the worse the angle distribution, except at λ = 1, where the optimal is consistently between h = 1 and h = 1.5.

C. Densities

By using N pT reference simulations and by computing the probability to generate not only the exact grain positions, but also the exact box dimensions, the STM method provides force field parameters that could lead to a good mass density. However it was observed in ref. 40, that the STM method fails to reproduce density at λ > 1. The interpretation is that the STM method has to compromise between reproducing the grain motion and the box size evolution and the force field model is not able to produce both quantities accurately at the same time. By using a more realistic, anisotropic model, we could expect that the reproduction of both quantities could be improved. The result is shown in Figure 6.

It is obtained, that at λ = 1, density is accurately reproduced whatever h. The coarse-grain density remains within 2 % of the MD target with a slight increase of the discrepancy at intermediate h.

At λ = 2, the density is quite sensitive to h and in-FIG. 6. Equilibrium density at T = 300 K and p = 1.0 bar. The agreement with the MD target is excellent at all h for λ = 1, within 3 % of the target for h ≤ 2.5 (optimal at h = 2) at λ = 2 and for h ≥ 3.25 (optimal at h = 4) at λ = 3.

creases from h = 1.25 on. It remains within 3 % of the MD target up to h = 2.5 and crosses the target close to h = 2. From this h-value on, the deviation from the MD target value is increasing with h.

At λ = 3, the evolution of the curve is completely different and the density is always underestimated by the coarse grain model. The agreement steadily improves with increasing values of h. Optimum values of h for reproducing the density within 3 % start from h = 3.25.

To complete this density-dependence on h, it should be noted that this property is also sensitive to the range of interactions and as a result on the cutoff of pair interactions. But, this effect is not presented in this work since the choice of the cutoff value was based on the probability P.

D. Compressibility

In order to compute the compressibility χ T , each model was used to produce trajectories at different pressures (p = 0.1, 0.5, 1.0, 5.0, 10.0 bar) and the equilibrium density ⟨ρ⟩ was obtained in each case. Compressibility was then obtained by fitting a straight line to the ln⟨ρ⟩ vs. p data points. Indeed,

χ T = ∂ ln⟨ρ⟩ ∂p T (5)
Example of such fits are shown in Figure 7 at λ = 3. The confidence interval is computed using Student's law with a confidence level of 95 %. The compressibility as a function of h at different λ is shown in Figure 9. As expected, compressibility increases with the degree of coarse graining λ. Unfortunately, overall, introducing anisotropy does not resolve the compressibility issue, i.e. an excessive density-dependence on pressure. The curves are quite noisy, but the trend shows that anisotropy only slightly improves compressibility at small h. This means that overall the interactions are not significantly less soft by using an anisotropic model. Interestingly, the slight improvement at small h correlates with the improvement of the model likelihood (see Figure 2).

The reason why the improvement in compressibility is not as great as expected may be due to the simplifying assumptions of the BOA model. Although realism is indeed improved by using anisotropic grains, as evidenced by Figure 2, the model could be even more realistic (at the expense of more computation time) by considering torques (see Figure 8), by using triaxial ellipsoids (or even dynamic shapes) and by letting the main axis of the grain depart from the strict backbone orientation.

We also note that despite the significant deviation with the MD reference by a factor 3, the CG model at λ = 1 remains fairly incompressible at all h because χ T is still very small.

E. Bond crossings at the mesoscale

Real polymer chains are not able to pass through each other due to the excluded volume. This leads to the formation of entanglements, chain confinement inside tubes, etc. In atomistic simulations, the chemical bonds are small and the pairwise repulsion is hard at short distance. This prevents bond crossing from occurring. In order to detect possible bond crossings during the CG simulations, each pair of bonds was tracked between each consecutive time step and the algorithm from ref. 59 was used. At λ = 1, the interactions are still hard enough and the bonds short enough to prevent any crossing. At least in the equilibrium simulations that we have run, no bond crossings were observed. At λ = 2, the interaction becomes softer and chain crossings are detected. However, by increasing h, the number of crossings drops by more than 90 % (see Figure 10). At h > 2, anisotropy also reduces chain crossing, but the effect is less spectacular. The reduction is by about 50 % at λ = 3, 30 % at λ = 4 and 15 % at λ = 5. FIG. 10. The relative number of bond crossings in terms of the anisotropy factor h, and at different coarse-graining levels λ. For isotropic grains (h = 1), the average number of bond crossings per bond and per nanosecond are 0, 0.069, 1.461, 3.213 and 3.818, for λ = 1, 2, 3, 4 and 5, respectively. Note that no bond crossings are detected at λ = 1, whatever h.

F. Diffusion

The STM method analyzes the short-time dynamics of the full system to optimize the CG force field. Now we focus on the long-time dynamics of the chains. In order to characterize the long-time dynamics, we computed the mean square displacement (MSD) as a function of time.

To do so, we simulated a system of 20 long entangled chains of 900 monomers each. As is often the case, we found that the MSDs of the monomers and of the whole chains were too fast at the DPD scale compared to their values at the MD scale. See Figure 11. This is partly explained by the possibility of unphysical bond crossing and tube escaping. But other causes must contribute, since the same effect is observed at λ = 1 but less pronounced than at higher levels of coarse-graining, where no bond crossings occur.

When h increases, the MSD slows down at all level of coarse-graining. At λ = 1, the MSD eventually becomes too slow beyond h = 2, while at λ = 3, the MSD always remains too fast. At λ = 2, where the reduction of bond crossings is the most significant, the agreement between DPD and MD MSDs is quite satisfactory at h = 4. Please note, that

• we have simulated a system, where each chain is entangled with its own periodic image, which is not ideal. We are however confident, that similar results would have been obtained if we had been able to simulate much bigger systems. • the DPD friction used is optimized independently for each anisotropy following the STM method and results in friction values equal within a few percent. Although MSD is very sensitive to this value, the important differences observed here in the MSD cannot be caused by such small differences in the friction.

• the optimal anisotropy for reproducing the shortterm dynamics (ie. 1 time step, 50 fs) is actually given by the likelihood shown in Figure 2. We observe that the optimal anisotropy is not the same when considering short-term or long-term dynamics.

IV. CONCLUSION

The development of coarse-grained models is a research activity in full development because it is difficult to design lower resolution models capable of rendering the chemical nature of the polymer.

We have shown that anisotropy can improve significantly the model realism. Depending on which property is the most important, different choices could be made. The optimal values of h are summarized in Figure 12. It appears that the optimal anisotropy is not unique for all criteria. For example, at λ = 2, the optimal dynamics (likelihood) is obtained with h = 1.5 while the optimal local structure (RDF) is obtained with h = 3.25. As a consequence, it is not possible to recommend a globally "best" anisotropy. Instead, the practical choice of the anisotropy should result from a trade-off between several criteria.

Of course, computation time is also a criterion of primary importance. Figure 13 shows how it increases with the anisotropy factor. The reason why it increases is that the number of neighbor particles interacting with anisotropic particles is greater because the distance r between particles may be greater than the effective distance R.

All things considered, in the particular case of cPB, we may recommend to use 1 ≤ h ≤ 1.75 at λ = 1, 2 ≤ h ≤ 2.5 at λ = 2 and h = 4 at λ = 3. The greatest benefit of using anisotropic grains would then be a better match of structure at λ = 1, much fewer bond crossings at λ = 2 and a much better density at λ = 3.

For going beyond the BOA model, we think that the compressibility issue may be solved in all systems by using the density-dependent potentials 37,38 . The STM method could be used to obtain the interaction parameters in this case. At λ = 1, the model is already quite good. Perhaps using non-spheroid ellipsoids could improve it even further. At high coarse-graining level, using deformable grains is probably the way to go. The bond crossing issue remains a serious problem beyond λ = 2. The STM method could also be used to parameterize segmental repulsions as in Refs 42,44.

V. SUPPLEMENTARY MATERIAL 1) Force field parameters at the MD scale. 2) STM likelihood at λ = 4, 5. 3) Functional form of the interactions at the CG scale. 4) Potentials at λ = 4, 5. 5) Bond and angle distributions with their corresponding errors. The used polymer systems contain 18 000 chemical monomers (18 000, 9000 and 6000 CG particles for λ = 1, 2 and 3, respectively).

FIG. 1 .

 1 FIG. 1. Mapping of one cPB chain, which contains 60 monomers of 1,4-cis-butadiene. From top to bottom, 1, 3 and 5 monomers are mapped into a spherical CG bead. The terminal beads are similar to the others, except that their molar mass is increased by 1 g mol -1 due to the extra hydrogen atom.

FIG. 2 .

 2 FIG. 2. Average values of the negative logarithm of the probability P, computed during the STM process. For each coarsegraining level λ ≤ 3 such parameters are evaluated at different anisotropy factors h, and at different cutoff radius rc for conservative interactions. The minimum in the curve indicates the most likely value of h, ie. the anisotropy factor with which the probability to generate the reference trajectories is the highest. Plot c also shows that the most likely cutoff at λ = 3 is rc = 20 Å. The corresponding parameters for λ > 3 are given in the Supporting Information.

FIG. 3 .

 3 FIG.3. Pairwise interaction potentials at different coarse-graining levels and different anisotropy factors. For the coarsegraining levels λ > 3, the corresponding potentials are shown in the Supporting Information. Note that the abscissa is the orientation-dependent effective distance R, not r.

FIG. 4 .

 4 FIG.4. Radial distribution functions. The structure is better matched using a given h > 1 at each coarse-graining level λ. See Figure5for a quantitative measure of the match.

FIG. 7 .

 7 FIG. 7. Fitting in straight line of the data points, obtained by the evolution of ln⟨ρ⟩ in terms of the pressure p.

FIG. 8 .

 8 FIG.8. Following the BOA model, the anisotropic grains exert non-central forces on each other (black arrows), but there is no direct torque tending to reorient the grains (white arrows).

FIG. 9 .

 9 FIG. 9. Isothermal compressibility χT at T = 300 K versus the anisotropy factor h. The red dotted line represents the reference value (6.1 × 10 -4 MPa -1 ) obtained in MD simulations.

FIG. 11 .

 11 FIG. 11. Mean Square Displacement (MSD) as a function of time for different anisotropy h and CG level λ. The system consists of 20 chains of 900 chemical monomers, simulated in the N pT ensemble at 300 K and 0.1 MPa. The bottom curves are the MSD of the whole chain, while the top curves are the MSD of individual monomers. Compared to MD values, MSD is always too fast at λ = 1 and slows down when h increases.

FIG. 12 .

 12 FIG.12. Optimal value of the anisotropy factor h depending on the level of coarse-graining λ and the chosen criterion. A single optimum value for all criteria cannot be given. Please note that the ranges may be subjective.

ACKNOWLEDGMENTS

This work was performed in SimatLab, a joint publicprivate laboratory dedicated to the modeling of the polymer materials. This laboratory is supported by Michelin, Clermont Auvergne University (UCA), CHU Clermont and CNRS. Computer simulations were performed using the facilities of the supercomputer Mésocentre Clermont Auvergne.