{'fraction_non_alphanumeric': 0.02891359692419598, 'fraction_numerical': 0.021535823764742557, 'mean_word_length': 4.503216583273767, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A coarse-grained model for capturing the helical behavior of isotactic polypropylene

Understanding the process–property relations of helical polymers using molecular simulations has been an attractive research field over the years. Specifically, isotactic polypropylene still remains a challenge for current computational experimentation, as it exhibits phenomena such as crystallization that emerge on large spatial and temporal scales. Coarse-graining is an efficient technique for approaching such phenomena, although previous coarse-grained models lack in preserving important atomistic and structural details. In this paper we develop a new coarse-grained model, based on the popular MARTINI force field, that is able to reproduce the helical behavior of isotactic polypropylene. To test the model, the predicted statistical and structural properties (characteristic ratio, density, entanglement molecular weight, solubility parameter in the melt) are compared with previous simulation results and available experimental data. For the development of the new coarse-grained force field, a single unperturbed chain Monte Carlo algorithm has been implemented: an efficient algorithm which samples conformations representative of a melt by simulating just a single chain.


Introduction
Polypropylene is a material that has been widely used for over 50 years. Since the first polymerization of propylene in 1951 by Paul J. Hogan and Robert L. Banks, 1 polypropylene has attracted attention due to its versatility and its complicated semi-crystalline morphology. Fibers, 2,3 films 4,5 and automotive peripherals 6,7 are some of its applications. In order to get insights into the processability, drawability and processproperty relations of the product during manufacturing, both melt and semi-crystalline polypropylene morphologies have been systematically studied. [8][9][10] Important insights into both structural and dynamical properties of polypropylene can be provided by modern molecular simulations, using molecular-dynamics (MD) [11][12][13] and Monte Carlo (MC) techniques. 14,15 Despite considerable progress, phenomena such as the development of polymorphism 16 are still not accessible via current simulation computer power. Furthermore, the understanding of complex and slow relaxations in polypropylene melts still remains a challenge.
A promising strategy for exploring in silico the large spatial and temporal scales of polypropylene structural and dynamic behavior is coarse-graining. In a coarse-grained model, a number of heavy atoms is combined into a single coarse-grained interaction site (bead). In several coarse-grained approaches used so far-bottom-up, [17][18][19][20] top-down [21][22][23] and hybrid 24,25 -two main characteristics can be distinguished: (a) mapping the atomistic structure onto the coarse-grained level and (b) development of a coarse-grained force field that describes the interactions among beads. Back in 1995, the first attempt at coarse-graining polypropylene was made by Schweizer et al., 26 who tried to implement PRISM theory 27 using the rotational isomeric state model (RIS) on different coarse-grained schemes. Then, Mattice and coworkers 15,28 developed a coarse-grained iPP model which was placed on a high coordination lattice, with each site corresponding to an iPP monomer. Later, based on the MARTINI force field, a coarse-grained scheme specifically for polypropylene was proposed by Panizon et al., 29 with each bead comprising a methyl group, a methine group and two half-methylene groups (one bead per repeat unit). Recent studies use even coarser schemes, choosing a Kuhn monomer to represent a bead. 30,31 The challenge in polypropylene coarse-graining is connected to its structure. Specifically, polypropylene is a branched vinyl polymer which exhibits a tendency to form helices. These helical conformations are drastically affected by the stereo-regularity of the branch. For instance, while isotactic polypropylene blocks can adopt a helical conformation upon cooling down below the melting point, atactic blocks remain amorphous even at low temperatures. Short helical configurations are also present in polypropylene melts. Upon coarse-graining, as the resolution becomes less fine, the methyl branches disappear and their influence on conformation must be accounted for indirectly. Until now there is no off-lattice coarse-grained model that describes polypropylene's helicity.
The purpose of the present study, which focuses on isotactic polypropylene (iPP), is twofold: first, we introduce and check the stability of the iPP helical structure based on existing coarse-grained representations. Secondly, we implement a new coarse-graining method and test the force field produced for iPP. An iPP helix can be represented at the coarse-grained level with distinct handedness and its structure can be preserved. To develop such an iPP coarse-grained model that preserves helicity, a hybrid approach is followed, in the sense that united-atom simulations are carried out for the parametrization of the bonded interactions, while the non-bonded interaction potential is tailored to reproduce experimental results. The novelty of the method lies in that a single unperturbed chain Monte Carlo algorithm is employed at the united-atom level as a reference in order to extract the bonded coarse-grained interactions. The new coarse-grained model has been developed having crystallization in mind. Simulations of crystallization will require a separate study and are not within the scope of this paper. Only some preliminary results are provided in the Discussion section to show that the model has much potential for providing insight into the crystallization of polypropylene.
The paper is organized as follows. In Section 2 we present the models and the simulation techniques that are used in the present study. In Section 3 we develop the new HELical Isotactic Polypropylene Potential (HELIPP) at the coarse-grained level, and validate it by comparing against other existing experimental and simulation data. Finally, in Section 4 we discuss the limitations of the new HELIPP and the possibility to tackle the remaining iPP modelling challenges.

Models
Isotactic polypropylene (iPP) has been simulated both at the united-atom (UA) and the coarse-grained levels (CG). The UA representation, in which carbons and hydrogens are fused into one single atom, as shown in Fig. 1a, is described by the TraPPE-UA force field that was initially developed by Martin and Siepmann 32 and then improved by Pütz et al. 33 by introducing an improper dihedral potential that preserves the correct stereochemistry. In the UA model, the non-bonded interactions are calculated for atoms separated by more than three bonds, as TraPPE-UA excludes 1-4 interactions. At the CG level, a single interaction site (bead) comprises a methyl, a methine and two half methylene groups that are shared between adjacent beads ( Fig. 1b). The center of the CG site is placed at the center of mass of these groups. The latter representation is invoked in the MARTINI force field 21 which was refined specifically for polypropylene by Panizon et al. 29 In this MARTINI model, the nonbonded interactions between the first and the second neighbors are excluded.
Monte Carlo and molecular-dynamics simulations were performed, as explained later, for a single iPP chain (both in a perfect helical structure and coiled), as well as for iPP melts. All the systems, except for the perfect iPP helical structures, were generated using the Amorphous Builder plugin of the Materials and Process Simulation (MAPS) platform developed by Scienomics. 34 The perfect iPP helical structures at the UA level were created based on diffraction data from Immirzi and Iannelli, 35 in the same way as it has been recently described by Romanos and Theodorou. 36 A perfect iPP helical structure, with three iPP monomers per helical turn, can be either left-handed (LH) or right-handed (RH). At the UA level a RH helix can be spotted as sequences of trans (1801) (CH 2 -CH-CH 2 -CH) and gauche À (À601) (CH-CH 2 -CH-CH 2 ) dihedral angles along a chain (Fig. 1c) resulting in an up arrangement, 37 while a LH helix is spotted as a sequence of gauche + (+601) (CH 2 -CH-CH 2 -CH) and trans (1801) (CH-CH 2 -CH-CH 2 ) dihedral angles resulting in a down arrangement. 37 For the dihedral angle definition, we make use of the IUPAC convention. 38 The two types of helices are enantiomers, meaning that one is the mirror image of the other. It is worth noting that high energy helices do exist and can be defined as a sequence of (trans, gauche + ) and (gauche À , trans). 39 In order to represent an iPP helical structure at the CG level, the UA RH and LH helices are mapped based on the mapping proposed by Panizon et al. 29 In Fig. 1d, it is clearly seen that the helical structure of the iPP chain persists at the CG level as well. The dihedral angle pattern for a CG helix should also be defined: a RH CG helix is denoted by a sequence of +1001 dihedral angles, and a LH CG helix is denoted by a sequence of À1001 dihedral angles. At the CG level there is no distinction between high energy helices and low energy helices.

Simulation methods
The first goal of the present study is to test if the existing MARTINI force field used for modeling iPP at the CG level 29 is able to capture the helical behavior of chains. For this purpose, the perfect iPP helical structures are simulated at both the UA and CG levels. The next step is to tailor the current CG force field aiming at the reproduction of melt structural properties both from a single chain MC simulation at the UA level and experimental measurements. The general strategy for developing the new CG force field is shown in Fig. 2. In Fig. 2b, a single UA iPP chain is presented, which will be simulated with MC by mimicking the interactions that a chain feels in an iPP melt (Fig. 2a). Mapping the produced sample of equilibrated UA conformations to the CG level (2c), CG bonded interactions can be extracted. CG non-bonded interactions are adjusted to experimental properties by using a geometric approximation, which is discussed in more detail in Appendix 1. Finally, the validation of the model has been carried out by performing both single chain MC simulations and melt molecular-dynamics simulations at the CG level.
2.2.1 Molecular-dynamics simulations. In order to check if helices of different handedness are energetically equally probable, single chain MD simulations have been performed at both UA and CG levels. In detail, single chain simulations of RH and LH helices were carried out at constant low temperature of T = 10 K using the velocity rescaling thermostat. 40 No boundary conditions were used. Coarse-grained helices were simulated under the same conditions. Then, for testing the ability of the existing MARTINI model 29 and the newly developed model to produce iPP melts of RH and LH helices in equal amounts, manychain iPP melts were simulated in the NPT statistical ensemble at T = 500 K and P = 1 bar using the Nosé-Hoover thermostat and the Parrinello-Rahman barostat. 41 In this case, periodic boundary conditions were used in all directions. The timesteps for the UA and the CG levels were 1 fs and 20 fs, 29 respectively. At the CG level, we used the same timestep as suggested in ref. 29, where the MARTINI implementation was discussed. Moreover, for the pair list generation the Verlet algorithm was used 42 with a cut-off distance of 1.4 nm. For the van der Waals interactions, a straight cut-off was chosen at r c = 1.2 nm at the UA level with tail corrections and at r c = 1.1 nm at the CG level with a shift of the potential to zero at the cut-off distance. All the simulations were performed using the Gromacs 5.0.7 43 simulation package. More details are given later in the Results section.
2.2.2 Single chain MC simulations. Single chain MC simulations have been carried out for a single iPP chain on both the UA and CG levels at T = 500 K. Four MC moves were implemented according to Tzounis et al., 44 each with a specific attempt probability. These MC moves are: (a) translational or single atom displacement (10%), 45 46 (c) pivot move (70%), 47,48 (d) flip branch move (10%). At the CG level, move (d) is omitted, as there is no branch, and the attempt probability of move (c) is increased by 10%.
On both levels (UA and CG), each force field describes explicitly all stretching, bending and torsional interactions. In a full melt simulation, non-bonded interactions would be calculated for atoms both along the same chain (intramolecular) and between different chains (intermolecular). Based on Flory's random coil hypothesis, 49 however, a chain in a melt behaves as if it is subject to local interactions only. For the purpose of sampling single unperturbed chain conformations representative of the melt, nonlocal intramolecular interactions between topologically distant segments and all intermolecular interactions can be omitted. Thus, to sample a single unperturbed chain, one needs to define precisely what constitutes local interactions. Tzounis et al. 44   Dn pair corresponds to the number of repeat units to the left and to the right, with which a specific repeat unit interacts. All the atoms belonging to the central repeat unit, marked with a rectangle, interact with all the atoms belonging to the neighboring repeat units, marked with ellipses. Dn pair = 3 was determined appropriate for polypropylene. 50 a way of defining local interactions by systematically varying the range of nonbonded interactions along a chain (Fig. 3). This range is set by a parameter Dn pair , which determines with how many neighboring monomers a specific monomer interacts and is estimated based on an empirical criterion. 44 Specifically, the value of Dn pair which maximizes the chain characteristic ratio is taken as representative of the unperturbed state. 44 By incorporating these local interactions we are able to sample conformations representative of a melt by simulating just a single chain. Dn pair was set equal to 3 for both UA and CG simulations based on previous work on polypropylene. 50

Testing CG MARTINI force field
In this section, the ability of the existing CG MARTINI force field of Panizon et al. 29 to produce iPP melts with equally probable RH and LH helices was tested. We started with single iPP helix simulations at the UA and CG level and we proceeded with iPP melt simulations at the CG level. UA RH and LH helices of 18 monomers each were simulated using MD at 10 K for 400 ps in total. The initial helical configurations were mirror images of each other, while a different set of initial velocities was assigned to each helical configuration following a Maxwell distribution. The reason the systems were simulated at such a low temperature is that it was not intended for any of the helices to lose their helical conformation. In particular, we aimed here at the computation of the dihedral energy of the two different helical types.
In Fig. 4a, the evolution of the dihedral energy of the RH and LH UA helices in time is plotted based on the united-atom representation. It is clearly seen that the dihedral energies for each helix fluctuate between the same energy levels. The fluctuations do not alter the helical structure. Next, the two helices were mapped onto a CG representation and a simulation of duration 400 ps was performed. In this case, the dihedral energy for the RH and LH CG helices came out significantly different (Fig. 4b).
Even though the helicity is preserved, despite fluctuations, the RH CG helix is energetically higher than the LH helix. Therefore, the MARTINI force field 29 favors one helical type (the LH helix) over the other.
To illustrate this further, we performed an MD simulation starting from an iPP CG melt of 40 chains of 50 monomers per chain. The iPP melt was first energy-minimized using the steepest descent algorithm 51 until the maximum force was smaller than 100 kJ mol À1 mm À1 . The algorithm converged after approximately 500 steps. Subsequently, a NVT simulation followed for 50 ns and, finally, a NPT simulation was carried out for 1000 ns. The final 500 ns were used for production. The simulation temperature (500 K) is well above the melting temperature for iPP oligomers that has been calculated at 411 K. 52 The simulated radius of gyration (1.48 nm) and density (1048 kg m À3 ) are in excellent agreement with the results of Panizon et al. 29 Note that this density deviates significantly from available experimental data (710-764 kg m À3 ). 53 In the CG iPP melt simulated with the MARTINI force field 29 we have analyzed the sequences of helical segments along a chain, in order to quantify the length of helical sequences formed spontaneously in the melt. More specifically, a CG segment was characterized as helical if |f + 1001| o 251 (LH) or |f À 1001| o 251 (RH), where f denotes the dihedral angle based on the IUPAC convention. 38 A sequence of such dihedrals along the chain is called helical sequence in what follows. In Fig. 5a, the distribution of lengths of helical sequences after 250 ns of NPT MD is plotted for the MARTINI 29 model. The percentages in parentheses indicate the fraction of the atoms belonging to either a LH or a RH helix. The initial configuration produced by MAPS had equally distributed helices. As the CG simulation proceeded, the population of the LH helices increased drastically above 40%, while that of RH ones decreased to about 15%. During the whole simulation, the helical sequences remained shorter than 11, the threshold beyond which ordering is initiated. 54 The overall behavior of the helices is also displayed in the probability density function (PDF) of different CG effective dihedral angles plotted in Fig. 5b. Two maxima can be seen, one at f = À1001, corresponding to the LH helix, and another one at f = +1001, corresponding to the RH helix. The f = À1001 peak is clearly higher and thus, by using MARTINI, 29 the LH helices in a melt are overpopulated.

New parametrization of CG dihedral and non-bonded interactions
The main goal now is to validate if the single chain TraPPE-UAbased MC simulation is able to reproduce correctly the structural melt properties and then use it as a reference for an improved CG model which preserves the equal population of both LH and RH helices. Error bars are included only if they are large comparable to the mean values. Otherwise, they are omitted. We simulated a UA single chain of 2000 iPP monomers using MC for 10 9 steps and the TraPPE-UA force-field, as was described in Section 2.2.2. To begin with, the characteristic ratio C N was calculated as follows: where hR 2 i is the mean-square end-to-end distance of the chain, n is the number of chemical backbone bonds and l 2 is the average squared chemical backbone bond length. C N was estimated at 6.4, while the experimental value reported by Fetters et al. 55 is 6.2. There is a slight overestimation, but the agreement is acceptable. Available melt iPP simulation results with a different force field reported C N = 5.49 AE 0.33. 56 The underestimation of C N in ref. 56 can be attributed to the short chain length of the system used there (152 skeletal bonds per chain). Following Brandrup et al., 57 the iPP chain radius of gyration R g was estimated as where C = 0.032 mm mol 1/2 g À1/2 and M is the molecular weight in g mol À1 . The resulting value for a 2000 units iPP chain, based on experimental data, 57 is 9.27 nm. The present simulation gives R g = 9.8 nm. Overall, the TraPPE-based C N and R g denote a slightly stiffer chain but overall the results are in a good agreement with experimental and theoretical predictions. In addition, the helical sequence length distribution for UA iPP is plotted (Fig. 6a), as also done by Yamamoto. 39 The iPP configuration was extracted after 5 Â 10 8 MC steps. It is important to note that the populations of both types of helices are equally distributed. Despite having a long chain of 2000 monomers, the helical sequences do not grow to more than 12 monomers per sequence. Thus, this UA model has no preference for any of the helical types, preserves the helical behavior, and reproduces the conformational properties of realistic iPP very well.
The next step was to use TraPPE-UA results as a reference in order to parametrize the dihedral effective potential for the new CG force field which would preserve the equal distribution of both LH and RH helices. For this purpose, the UA trajectory was mapped onto the CG level. From this mapped CG trajectory the PDF P(f) of dihedral angles was extracted and, subsequently, was converted to the dihedral effective potential, U(f) = ÀRT ln(cP(f)). The constant c was selected arbitrarily equal to 10001, as it does not have any effect on the simulation. Then, a Ryckaert-Bellemans function (eqn (5)) was fitted to U(f), as shown in Fig. 6b. In the same figure, the MARTINI dihedral effective potential U(f) is also plotted, where at f = À1001, corresponding to LH helix, the dihedral energy is significantly lower than the dihedral energy at f = 1001, corresponding to RH. The new CG force field comprises effective stretching, bending, torsional and non-bonded interactions, with the parameters presented in Table 1. The reader should note that any coupling between successive effective torsion angles is ignored in the CG model. In other words, the total CG effective torsional potential of a chain is modeled as a sum of terms, each depending on a single effective torsion angle.
Stretching interactions. The CG stretching interactions are described as: where r ij denotes the distance between two CG beads connected by a bond with force constant k b . The mean CG bond length l o was computed equal to 0.293 nm. The mean bond length used by Panizon et al. 29 was 0.298 nm, and we make use of the same k b and l 0 as they did. 29 Bending interactions. The CG bending interactions are described by a Restricted Bending (ReB) potential as proposed by Bulacu et al., 58 where y ijk is the angle that is formed by three CG beads i, j, k, and k y is the bending constant. The value of y o is 1191, the same as in the iPP MARTINI model. 58 The reason for selecting the specific potential is the advantage provided by ReB in using longer timesteps and avoiding numerical instabilities that are often encountered in very flexible coarse-grained molecules. 58 Torsional interactions. The CG effective dihedrals were parameterized based on the single chain Monte Carlo Simulations explained previously in Section 2.2.2. For the fitting, we choose a Ryckaert-Bellemans function, n¼0 a n cos f ijkl n ; ( where f ijkl is the dihedral angle formed by four CG beads i, j, k, l. The parameters extracted from the fitting to the dihedral effective potential curve are noted by a n . Non-bonded interactions. For the representation of the nonbonded interactions in the CG model, a Lennard-Jones potential was used, where r ij is the distance between two interacting particles i and j, s is the distance at which the energy is zero and e is the depth of the potential well. We use the same e as Panizon et al., 29 while we set s = 0.53 nm, based on a geometric approach (see Appendix 1), in contrast to s = 0.43 nm used in MARTINI. The non-bonded interactions are calculated both intermolecularly and intramolecularly for CG beads separated by at least four CG beads (monomers).  (3)) l o = 0.298 nm, k b = 48 000 kJ mol À1 nm À2 Bending interactions (eqn (4)) y 0 = 1191 k y = 78 kJ mol À1 rad À2 Torsional interactions (eqn (5)) a 0 = À8.08655 kJ mol À1 , a 1 = À2.88516 kJ mol À1 , a 2 = 16.61891 kJ mol À1 a 3 = 0.79249 kJ mol À1 , a 4 = À8.79016 kJ mol À1 , a 5 = 1.33047 kJ mol À1 Non-bonded interactions (eqn (6)) e = 2.625 kJ mol À1 , s = 0.53 nm

Validation of the HELical isotactic polypropylene potential (HELIPP)
To test the HELIPP CG force field, a single unperturbed CG chain of 2000 chemical monomers was simulated by Monte Carlo for 10 9 steps. The optimal value for Dn pair was estimated equal to 3. The characteristic ratio C n vs. the number of real chemical bonds N b is plotted in Fig. 7. It is obvious that C n has clearly reached a plateau value for N b 4 3000. The plateau value of C n calculated by HELIPP equals 5.7, slightly lower than the value measured using TraPPE-UA (6.4). By performing single chain MC simulation using the existing MARTINI force field, 29 the estimate was 8.9. So, HELIPP reproduces better the stiffness of the iPP chain than the existing MARTINI CG model 29 and yields a characteristic ratio comparable to the experimental results. 55 All the simulations of the characteristic ratio were carried out at 500 K, while the reference temperature for the experimental measurements was 460 K. 9 The radius of gyration from single unperturbed chain CG simulations equals 9.23 nm, and is also close to the theoretical calculations. 57 Next, the PDF of CG effective valence angles is plotted in Fig. 8a. Looking at the HELIPP curve, a peak is observed at 1191, while the TraPPE-UA curve has a maximum at 1091, which is the most probable one, and a second weaker shoulder at 1451. As the ReB potential used for describing bending interactions 58 is not able to reproduce a double-peak curve, the two curves do not match each other. The HELIPP curve actually replaces the UA bending angle distribution with a smoothed, unimodal approximate distribution.
In Fig. 8b, the PDF of CG effective dihedral angles is also plotted. In both curves, two maxima can be distinguished. Each peak represents a helical type, both are equally probable. In the regions [À1801, À1501] and [1801, 1501], representing the trans states, the HELIPP curve exhibits a slightly higher probability. On the contrary, in the region [À301, 301], which refers to the cis state, the HELIPP curve indicates a lower probability. CG beads, which have a parameter s from the LJ interactions larger on average than the same parameter of the united atoms they are constituted from (see Appendix 1), experience slightly larger repulsion when approaching each other, which leads to the slightly less probable cis state.
Furthermore, we performed MD simulations of a CG melt of 40 chains of 500 beads each using this new HELIPP force field. In order to equilibrate the system, an energy minimization was undertaken until the maximum force was smaller than 10 kJ mol À1 nm À1 , followed by a NVT MD of 5 ns and a NPT MD of duration 5 ms. The simulation temperature T = 500 K was well above the melting temperature of high molar mass iPP (T = 460 K) measured by differential scanning calorimetry. 9 The simulation pressure P was set to 1 bar using the Parrinello-Rahman barostat. 41 The energy minimization converged after 25 000 steps.
We compare the MD simulation to the Single Chain Monte Carlo simulation with the same HELIPP force field. The PDF for  CG effective valence and dihedral angles extracted from the MD simulation match perfectly (not shown here) the ones derived from CG single unperturbed chain MC sampling and shown in Fig. 8. The characteristic ratio from the CG melt MD simulation was estimated equal to 5.87 AE 0.09, close to the value 5.7 from the single CG chain MC simulation. Moreover, the density was measured as 730 kg m À3 , in a good agreement with experimental results (710-764 kg m À3 ). 53 The helical sequence length was also studied. Both types of helices are now equally probable with approximately 30% of the beads belonging to RH helix and 30% to a LH helix during the whole simulation. Neither of them grows longer than 12 monomers per helical sequence. So, HELIPP is able to reproduce LH and RH helices in equal amounts and that holds even for rather long simulations on the ms scale.
Subsequently, using this new coarse-grained model, the entanglement molecular weight M e was calculated by performing a topological analysis on the trajectory produced during the last 1 ms of the simulation. For the topological analysis, the CReTA (Contour Reduction Topological Analysis) algorithm 59 was used. This algorithm maps an entangled polymer system to its corresponding system of primitive paths, having the same topological constraints as the original system. M e was computed equal to 4432 AE 121 g mol À1 , whereas in experiments carried out for iPP melt M e was estimated between 5100-5500 g mol À1 60 and 6900 g mol À1 . 61 So HELIPP produces slightly more entangled system in comparison to the experiments, but no comparison was made with other simulations, as no data were available.
Finally, we have calculated the cohesive energy of the iPP melt which can be expressed by the Hildebrand solubility parameter d. 62 d can be obtained as where E coh is the cohesive energy, n is the number of molecules, V cell the volume of the simulation box, N A Avogadro's number, U i,LJ the intramolecular LJ interactions energy of a single chain and U cell,LJ the intramolecular and intermolecular LJ interactions energy in the simulation box. The angular brackets denote ensemble averaging. In order to calculate d for different temperatures, we performed multiple NPT MD simulations for 10 ns for a system of 40 chains of 50 monomers per chain using the HELIPP and TraPPE-UA force fields. The simulation temperatures were T = 400, 450 and 500 K, close to the temperature at which the HELIPP force field has been developed. The temperature dependence of the solubility parameter is plotted in Fig. 9

Discussion
In this study we have developed a new CG force field for iPP melts, named HELIPP, which reproduces structural properties sufficiently well and, most importantly, is able to capture the different handedness of the iPP helices. HELIPP is remarkably simple and quite friendly to be adopted in simulations, as only one type of torsional angle is used and non-bonded interactions are described by a single LJ s and a single LJ e value. This force field could be possibly extended to other ordered polymer structures, as the approach we have developed is applicable not only for iPP melts. Furthermore, we have presented a method for developing such a CG model based on a single unperturbed chain MC algorithm and we discuss the possibility of using this method for polymers of rather arbitrary structure.
The representation of the RH and LH helices at the CG level and preservation of the handedness have been proved essential for the modeling of iPP melts. The existing MARTINI CG force field for iPP 29 is not able to reproduce correctly the two different helical types and overestimates the stiffness of the chains and the density. By introducing the CG iPP helices and tailoring the dihedral effective potential so that none of the helices is favored over the other, we managed to have a more accurate estimation for the iPP chain characteristic ratio, C N = 5.7. Moreover, by re-adjusting the Lennard-Jones s based on a geometric approximation, as discussed in Appendix 1, the density r = 730 kg m À3 was also improved. This force field was applied to an iPP melt but seems promising also for iPP semi-crystalline morphology. Below we present a snapshot of a lamella formation after an iPP coarse-grained system was submitted to stretch-induced crystallization using the HELIPP. Specifically, a system consisting of 20 chains of 1000 CG monomers each was equilibrated in the NPT statistical ensemble at T = 500 K and P = 1 bar using the Nosé-Hoover thermostat and the Parrinello-Rahman barostat. 41 Then, it was cooled down to T = 380 K with cooling rate 10 K ns À1 and stretched at the same temperature to 5 times its initial length in the N : e xx P yy P zz T statistical ensemble with stretching rate : e xx = 10 7 s À1 . Finally, it was annealed in the NL xx P yy P zz T statistical ensemble at T = 380 K for 400 ns (Fig. 10).
iPP crystals consist of lamellae which are formed by RH and LH helices in equal amounts. This characteristic is possible to maintain now at the CG level, although attention should be paid to two main aspects: (1) the presented force field has not been tested thoroughly for temperatures below the melting point. (2) At these low temperatures the reproduction of the different iPP crystal phases is crucial, and also has not been tested yet.
Single unperturbed chain MC simulations have been used for the fast and effective estimation of the stiffness of polymer melts but, as far as we know, not for the development of a force field up until now. The CG stretching, bending and torsional interactions can be accurately derived from a single chain MC simulation. In the present study we tailored only dihedral and non-bonded interactions, based on the approach implemented by Panizon et al. 29 For bending interactions we adopted the Restricted Bending (ReB) potential. In addition to the work presented above, we have performed CG simulations with tabulated potentials both for bending and torsional interactions which perfectly match the PDF of the CG effective valence and dihedral angles derived from TraPPE-UA single chain MC simulation (Fig. 8a). The results using tabulated potentials can be found in the ESI. † The main reasons we did not use these potentials for the production simulations were: (1) tabulated potentials can not be easily implemented in different simulation packages. (2) The maximum timestep we had to use with the tabulated potentials was 5 fs, while with ReB we use a timestep of 20 fs. For these reasons, we incorporated ReB and Ryckaert-Bellemans potential into HELIPP. As far as the nonbonded interactions are concerned, here we used a geometric approach which proved to be effective (see Appendix 1). This geometric approach is recommended to be tested also for other polymer systems. Another method that could be used for mapping non-bonded interactions is iterative Boltzmann inversion, 65 but, overall, the single chain Monte Carlo method can serve as a new coarse-graining method for mapping a variety of polymer melt models to the coarse-grained level.
Finally, we demonstrate the computational significance of the coarse-graining by comparing a UA and a CG simulations in terms of time. Coarse-grained models are far more efficient computationally than models with higher resolution. To check this, we simulated 20 iPP chains of 500 monomers per chain at the UA level using TraPPE-UA and at the CG level using the HELIPP. The simulations were carried out on a cluster made of 504 nodes comprised of 128 AMD Rome 7H12 processors, each with 2.6 GHz and 2 GB per core. The estimation made on 32 cores on 1 node for the CG performance was 4000 ns day À1 , while for the UA the performance was 100 ns day À1 . So, the CG iPP model performs 40 times faster than the UA iPP model.

Conclusions
The main scope of this research was to integrate the helical behavior of iPP into a coarse-grained model. For this reason, we developed a model that accounts indirectly for the influence of the iPP branch which is translated to different helical conformations along the chains. We showed that, even in melts, where the helical segments are short, the incorporation of helicity improves significantly the structural properties of iPP. In order to develop a new CG force field (HELIPP) that preserves helices, we implemented a coarse-grained method based on a single chain Monte Carlo algorithm which proved to be effective. Using the new HELIPP force field, long iPP chains were sufficiently relaxed in the melt by ms-long simulations and the characteristic ratio was accurately measured. In this way, slow relaxation above the melting temperature has been circumvented. In the future it would be interesting to extract dynamical properties as well from the new CG model. Furthermore, given that HELIPP captures the helical behavior of iPP which is essential for structural ordering and performs really well in terms of CPU time, studies of crystallization by molecular simulations using the new HELIPP model is an immediate objective for future work.

Conflicts of interest
There are no conflicts of interest to declare.

Appendix 1
Calculation of the coarse-grained Lennard-Jones r CG parameter In order to calculate the radius s CG at which the non-bonded LJ potential energy for coarse-grained particles is zero, we followed a geometric approximation. Given a bunch of united- Fig. 10 A lamella snapshot of a CG system of 20 chains of 1000 CG monomers each using the HELIPP after a stretch-induced crystallization simulation.
atom LJ spheres with specific radius s, the problem is to find the radius s CG of a coarse-grained LJ sphere with its center lying at the center of mass of these united atoms, that is approximately tangential to the sphere of each individual united-atom. In the present case, the united atoms that form a coarsegrained bead are: a methyl group (CH 3 ), a methine group (CH) and two half-methylene groups (CH 2 ). These united atoms can be placed in a Cartesian coordinate system with the methine group lying at the origin. The positions of the remaining united atoms are calculated according to the bond length and valence angles given by Pütz  All coordinates have been calculated in nm. In Fig. 11, the radius of the coarse-grained LJ sphere s CG/CH2 that is tangential to the LJ sphere of the methylene group and has its center at the position of the coarse-grained bead is calculated, s CG/CH 2 = 8r CG À r CH 2 8 + s CH 2 = 0.56 nm, where r CG is the position of the coarse-grained bead (black dot in Fig. 11), r CH2 is the position of the methylene group (green dot in Fig. 11) and s CH2 is the distance at which the Lennard-Jones potential energy for the methylene united-atom is zero (green arrow in Fig. 11).
Similarly, s CG/CH 3 = 0.49 nm and s CG/CH = 0.51 nm. All the distances at which the LJ potential becomes zero for the united atoms are given by Pütz. 33 Finally, s CG is calculated as the average of: ¼ 0:53 nm: Fig. 11 Lennard-Jones spheres of the united atoms and the coarsegrained bead. The coarse-grained LJ sphere with radius s CG is sketched in grey color, the methyl united-atom LJ sphere with radius s CH 3 in blue color, the methylene united-atom LJ sphere with radius s CH 2 in green color and the methine united-atom LJ sphere with radius s CH in red color. Each colored dot represents the center position of either the united atoms or the coarse-grained bead. (a) Profile view (b) side view.