Patrice Malfreyt

Keywords:

This doctoral work

Chapter 1 Introduction

Polymers are used in various industrial applications like plastic, textile, medical, food and beverages, cookware, automobile, aerospace, and other engineering fields. Considering an industrial world, the usage of polymers is of great importance. For instance, in the automobile industry polymers are used as part of advanced materials to reduce the overall weight of a vehicle. Another instance is the tire industry, where natural and synthetic rubbers are used to manufacture tires. Looking at the competitive edge like increased wear resistance, slower aging, reducing fuel consumption, increase tire life, and so on, Michelin along with many other tire companies such as Bridgestone, Goodyear, Continental invest in polymer research to meet the rising demands of the customer needs. Therefore, material designers are motivated in developing high-performance materials making it interesting to pursue this topic as research.

Polymer Structure

Polymers are long chain molecules that are composed of repeating units called monomers. Some polymers are found in nature, such as cotton, proteins, cellulose, silk, and rubber, while others are man-made or synthetic polymers, like polystyrene, polyethylene, and nylon. For instance, the primary constituent of natural rubber is polyisoprene, which can also be produced synthetically. Synthetic polymers are created via polymerization, which is the process of combining smaller molecules into a long chain held together by covalent bonds. A polymer can be composed of single or multiple types of repeating units. If a polymer is made of a single type of monomer, it is known as a homopolymer, e.g., polystyrene, which is composed of styrene units. When a polymer is made of more than one species or type of monomers, it is called a copolymer. The polymerization of different types of monomers is referred to as copolymersization. To discuss the different structures of copolymers more easily, we consider two monomers labeled by A and B:

• Block copolymers are made of repeating units in blocks like ....-A-A-A-B-B-B-.... . An example is a poly(styrene-butadiene-styrene) or SBS rubber, which is composed of a long chain of polystyrene, a long chain of polybutadiene in the center, and another long chain of polystyrene.

• When monomers are arranged in a pattern like A-B-A-B-A-B-, they are called alter-nating copolymers. Styrene maleic anhydride is a synthetic polymer made of styrene and maleic anhydride monomers.

• If the distribution of monomers follows some statistical law, they are called statistical copolymers. Random copolymers are a special case of statistical copolymers, where the probability of finding a given monomer at any given site in the chain is independent of the nature of the neighboring units.

• Graft copolymers are made of one or more side chains connecting to a main chain. The composition of the main and the side chains can be different, like polystyrene grafted onto polybutadiene.

Both homopolymers and copolymers can be classified into linear, branched, and network systems. Each type of these polymeric classification results in unique properties (e.g., mechanical, thermal, rheological, electrical).

A linear homopolymer is a chain with sequentially repeating monomers. An atomic model of a linear polymer is shown in Fig. 1.1, where one monomer is comprised of 13 atoms and the end monomers have 14 atoms, which includes an extra hydrogen atom. Figure 1.1: An atomic representation of a linear polymer with 80 monomers, where one monomer has 13 atoms of CH 2 =C(CH 3 )-CH=CH 2 and the end monomers have an additional hydrogen atom. The hydrogen atoms, carbon atoms with three, carbon atoms with two, carbon atoms with one and carbon atoms with no attached hydrogen atoms are depicted in blue, red, pink, yellow and green colors, respectively.

The short-range forces like bond and angle interactions act between the bonded atoms, whereas the long-range force like the van der Waals interactions act between non-bonded atoms. Examples of linear polymeric systems include polyesters, polyvinyl chloride and polyvinyl alcohol. Block and random copolymers can also be categorized as linear polymers, except that they are made up of different species of repeating units.

In molten state, the polymers can be molded and remolded by heating in industrial processing methods like injection molding, compression molding, calendering, and extrusion [START_REF] Macdermott | Selecting thermoplastics for engineering applications[END_REF]. Upon fast cooling the melt may form a glass or a glassy-state, which means that the polymers remain in an amorphous-state (random-coil-like) as in the melt, but the mobility of monomers is significantly reduced and the system becomes rigid [START_REF] Mattsson | Fluids, Colloids and Soft Materials: An Introduction to Soft Matter Physics[END_REF]. This generally occurs below a temperature called glass-transition or glass-liquid temperature (T g ). It is illustrated in Fig. 1.2 as the temperature at which the slope of the cooling curve changes. A decrease in temperature causes the particles to have less kinetic energy and reduced mobility, which results in particles that are close to remain even closer to each other and hence the volume decreases while cooling. When the polymeric system is cooled at a relatively slow rate, it may form a semi-crystalline phase, which has partially ordered regions surrounded by amorphous regions. The temperature (or range) at which the system transforms into a semi-crystal from the melt, is called the crystallization temperature (T c ) and it is depicted as the temperature at which the volume-temperature curve "drops" suddenly in Fig. 1.2. Generally, linear polymers are classified as thermoplastics, which means that they soften when heated, i.e., the intermolecular forces get weaker with increasing temperature and result in a melt-state with tangled chains. Branched polymers have single or multiple side branches protruding from the main chain forming a variety of structures, such as star-like, comb-like and dendrimer. The length of side branches may differ, giving rise to short and long-branched polymers. Branched polymers show a reduced efficiency in packing and generally result in low density systems [START_REF] Shrivastava | Introduction to Plastics Engineering[END_REF].

Volume / monomer

An example of branched polymers is low-density polyethylene (LDPE). An atomic representation of short branched LDPE is shown in Fig. 1.3. Similar to linear polymers, branched polymers are classified as thermoplastics and constitute nearly 76% of the global plastic usage [START_REF] Biron | Material Selection for Thermoplastic Parts[END_REF]. A cross-link is a connection between two polymers by covalent bonds. Several crosslinks can be formed, connecting different chains randomly to form a network. The covalent bonds in a cross-linked system are much stronger than van der Waals forces, leading to the increased stiffness and rigidity of the polymer system [START_REF] Shrivastava | Introduction to Plastics Engineering[END_REF]. A cross-linked system can be formed either with linear or branched polymers. An example of such a system is vulcanized natural rubber, i.e., when natural rubber is heated with sulfur, some double bonds between carbon-carbon react with sulfur, which allows to form a link or a bridge between different chains as shown in Fig. 1.4. Depending on the length and density of cross-links, the physical properties of the vulcanized rubber change [6]. A higher density of cross-linking can make a system hard and brittle, e.g., resin or plastic. However, the lightly cross-linked system has an increased flexibility and mobility like elastomers.

Elastomers are relatively long and entangled polymer network system. The schematic representation of a elastomer network with entanglements and cross-links are shown in Fig. 1.5(a) where Fig. 1.5(b) represents only entanglements and Fig. 1.5(c) depicts only crosslinks. A simplistic way to understand entanglements is to consider long polymers, that loop over themselves and other chains, forming a net-like structures.

The relatively low density of cross-links in the network system restricts the movement of chains when subjected to loading, thereby allowing the network to stretch to several times its original length [START_REF] Fried | Polymer science and technology[END_REF]. Elastomers are amorphous above T g , whereas below T g they form a glassy or semi-crystalline state depending on the rate of cooling. The mobility of polymers/segments of the polymers is more at higher temperatures leading to elastic behavior [START_REF] Mckeen | The Effect of Long Term Thermal Exposure on Plastics and Elastomers[END_REF]. Elastomers can be either thermoset or thermoplastic. However most elastomers are (S) n thermosets [START_REF] Mark | Applied Plastics Engineering Handbook[END_REF], which means these materials usually require vulcanization and they do not melt upon heating, for instance vulcanized natural rubber which is often used in tires and tubes manufacturing. Thermoplastic elastomers melt when heated because the cross-links formed are significantly weak, these materials are used in wire and cable insulation.

The mechanical properties of polymers are temperature-dependent. Fig. 1.6 shows a schematic representation of Young's modulus-temperature curve for a polymer system. At low temperatures, the inter-and intra-polymer movement is limited due to low thermal energy. As the temperature increases, the movement within and between the chains becomes higher. Based on the onset of motion on the scale of molecular displacements, the glass transition temperature can be determined [START_REF] Rudin | [END_REF]. Some of the common methods used to determine T g are Differential Scanning Calorimetry, Dynamic Mechanical Analysis, and Thermal Mechanical Analysis [START_REF] Zhang | [END_REF]. In the glassy region, polymers, in general, have high modulus and are brittle having the characteristic of glass. As the temperature increases, the modulus decreases. During the glass-to-rubber transition Young's modulus reduces by a factor of about 10 3 over a temperature range of 10 to 20 • C.

In the temperature interval of the rubbery plateau, the intra-chain movement is more but the movement of the inter-chain is still restricted. Such restriction can be due to the presence of cross-links and entanglements.

At even higher temperatures, the polymers can easily move. In the case of uncross-linked and thermoplastic polymers, the entanglements will slip and the system becomes a melt. In the case of elastomers, the cross-links present to restrict the flow of chains and the rubbery plateau persists.

Polymer crystallization

Crystallization can only occur in certain types of linear homopolymers, block copolymers, and elastomers. Such organization could be due to a reduction in the temperature of the system or due to the application of strain. When amorphous material like elastomers is subjected to stretching, the material can undergo a phase change and this phenomenon is called strain-induced crystallization.

Strain-induced crystallization (SIC) of polymers plays an important role in a variety of industrial applications, such as tires, fibers, and containers. The semi-crystalline structure induced upon straining is believed to be the origin of the high tensile strength and crack growth resistance of many polymeric materials like natural rubber, polychloroprene, polyethylene, poly(ethylene terephthalate), or polyvinyl alcohol (PVA) [12,13].

In elastomer networks and thermoplastics, crystallization plays a leading role during the processing stages thereby, giving the desirable mechanical properties of various industrial products such as tire tread bands for heavy-duty trucks and fibers or textiles used in bulletresistant vests and helmets. The structural properties and their relationship associated with the SIC of materials is a pretty old and well investigated topic. The number of studies carried out using different techniques is a sign that various topics are still a matter of debate. For example, it is unclear how the crystalline domains are changed upon deformation at the scale of crystallites [14,15], whether the evolution of the stretching ratio reaches a plateau during SIC [16], or what the impact is of the chemical nature of the polymer with regard to its kinetic for SIC.

The number of crystallites and the adjoined amorphous regions affect the mechanical, thermal, and optical properties. It has been hypothesized by Flory [17] that amorphous material between the crystalline regions is the reason for the high load-bearing capacity in elastomers during stretching. Since more than two-thirds of the global polymeric applications are covered by semi-crystalline polymers, this topic has received tremendous interest in the field of polymer research [18].

A common description of the crystallization process of polymers is based on the formation of crystals in simple particles, which includes the germination of nuclei referred to as primary nucleation and their growth into crystals which is termed as secondary nucleation [19].

Experimental observations

The study of polymer crystallization started in the early 19th century, when Gough and Page examined the elasticity of natural rubber [20]. However, they were not aware of the concept of stretch-induced crystallization while observing the elastic behavior in the natural rubber [18]. The phenomena of crystallization under strain was discovered by Katz in 1925 for a uni-axially stretched natural rubber using X-ray diffraction. Bekkedahl and Wood in the 1940s studied thermal-induced crystallization of natural rubber and they established that unstretched rubber crystallizes at temperatures between 20 • C and -50 • C, however at -25 • C the maximal crystallization rate was observed [21]. To examine the molecular structure during crystallization different techniques are used like transmission electron microscopy and wide-angle X-ray scattering [22].

It is worth noting that during the early 20th century Hermann Staudinger proposed that polymers are macromolecules and the atoms in these long chains are held by covalent bonds. Later, Herman Mark confirmed the hypothesis of long chains consisting of repeat units using X-ray crystallography technique. In 1953, Staudinger received the prestigious Nobel price in Chemistry for his contribution in the field of polymer science. The long chain systems cannot form fully ordered structures because of entanglements. The ordered regions of a semi-crystalline system are orders of magnitude smaller than the fully extended chain length. This observation led to the discovery of chain-folding phenomena using scattering experiments [23,24]. The folded chains stack together to form lamella, and later these lamellae grow radially into spherulites. Spherulites are not completely crystalline, but contain amorphous regions that connect different lamellae, as shown schematically in Fig. 1.7. The crystallization involves multi-scale process ranging from chain folding to formation of spherulites with the associated length-scales from Angstroms to micrometers and time-scales from femtoseconds to milliseconds [START_REF] Gedde | Fundamental Polymer Science[END_REF]. Rubber under different conditions, like uni-axial stretching and cyclic loading, exhibits SIC [22,[START_REF] Che | [END_REF]. A simplistic description for crystallization of rubber is that, upon straining, the entanglements along with the cross-links act like pivot points to align the chains and induce SIC. The different crystallites or lamellae are separated by amorphous sections of the chains. At large strains, polymers adjust the concentrated force in order to avoid the failure [28]. The presence of amorphous sections which is reducing between the oriented structures during straining is believed to be the key for high tensile strength and crack growth resistance in elastomers.

On the other hand, flow-induced crystallization is a phase transition behavior of polymer solutions or melts when exposed to an external flow like stirring or extrusion. Under flow the polymer solutions or melts yield an oriented morphology, where the chains stretch along the direction of flow and crystallize to form a fibrous core.

Spherulites can be formed from the stretched fibrous core or the stretched network, resulting in superstructures labeled as shish-kebabs or row-nucleated structures [START_REF] Gedde | Fundamental Polymer Science[END_REF]. A micrograph taken from the field emission scanning electron microscope of ultra-high-molecularweight polyethylene shows the formation of the shish-kebab morphology under shear in Fig. 1.8. Figure 1.8: Experimentally seen shish-kebab morphology of ultra-high-molecular-weight polyethylene upon shearing using a field emission scanning electron microscope. These shish-kebabs have been reported to have different lengths, which are of the order of µm and the outer diameter is approximated to be about 300 to 400nm [29] (Image from Ref. [29]).

Theoretical and Simulation studies

The nature of chain-folding has been described in several theories, here we discuss few of these models. Keller introduced an adjacent re-entry model for solution-grown crystals, wherein the chains fold back into the lamellae in a regularly ordered fashion [24]. It was experimentally observed that the crystals obtained from melt usually presented high mechanical strength and toughness. This observation led to the switchboard model by Flory [30], in which the chains do not re-enter the lamellae by regular folding, but rather reenter randomly and the lamellae must be connected to each other by amorphous sections of the chains, which are referred to as tie-chains [31]. The experimental observations indicated that the adjacent re-entry model is best suited for solution-grown crystals and the switchboard model describes the melt-formed crystals better. To capture the chain-folding phenomena in both types of crystal formation mentioned above, Hu et al. [32] performed Monte-carlo simulations of a single polymer for different chain lengths and computed the free energy. One of the main observations was that at different sections of the polymer, straight segments can be formed in a long chain simultaneously. This can be due to the constraints imposed by entanglements or the restrictions impinged by crystallites. In the intra-molecular nucleation model [31], the polymers may pass through multiple lamellae (tie-chains) or chain folding can occur in the same lamellae many times. The above-mentioned models are summarized in Fig. 1.9.

The widely accepted crystallization model is Hoffman-Lauritzen (HL) theory, which was presented as a kinetic theory. According to HL-theory:

• The crystallization starts with the formation of a stem (the straight sequences of chain which traverse the crystal), which is referred as a secondary nucleus

• After which growth takes place by adding stems in the vicinity of the first stem in a lateral fashion similar to the adjacent re-entry model and forming a layer of stems Repeating these steps leads to the growth of the crystal [18]. However, it is important to mention that each secondary nucleation can result in the formation of a new layer. The initial stage of crystallization referred as primary nucleation, which is associated with the germination of small nuclei before the formation of stems is not considered in this theory. The HL-theory provides an expression for the rate at which crystallites grow radially as a function of the degree of supercooling [START_REF] Gedde | Fundamental Polymer Science[END_REF]. Based on the growth rate of layers (g), which is the rate at which stems deposit next to each other laterally, and the rate (i) at which secondary nuclei are formed, HL-theory introduces three crystallites growth regimes [33]:

• Regime I: i g, usually at low degree of supercooling, where each layer or monolayer is completed before the secondary nucleation event occurs • Regime III: i g, high supercooling, where growth occurs as a result of many secondary nucleation events

The HL-model has been successful in providing a general guideline for crystallization kinetics. It predicts the crystallite growth rate and the crystal thickness [18,[START_REF] Gedde | Fundamental Polymer Science[END_REF]. The growth rate of crystals predicted by HL-theory is given by

G r ∝ exp - K g T c ∆T , (1.1)
where G r is the radial growth rate of spherulites and K g is a regime dependent kinetic constant. ∆T = T 0 m -T c is the degree of supercooling, where T 0 m is the equilibrium melting temperature and T c is the crystallization temperature. A detailed derivation of expressions of G r and K g are presented in Chapter 8 of [START_REF] Gedde | Fundamental Polymer Science[END_REF].

The crystal thickness L * c is given as,

L * c = C ∆T + δL c , (1.2)
where C and δL c are regime dependent constants. It is reported that crystal thickening is not considered in HL theory [START_REF] Gedde | Fundamental Polymer Science[END_REF].

Although the HL-model was initially developed for crystallization from solutions, it has been refined to be applicable to crystallization in melts [33]. In particular, the melting temperature T m of a semi-crystalline system is expected to vary linearly with the inverse lamellar thickness [START_REF] Muthukumar | Advances in Chemical Physics[END_REF] and this behavior is observed in many experiments. However, some experimental results show that HL-theory fails in predicting some morphological features such as [START_REF] Gedde | Fundamental Polymer Science[END_REF]:

• The assumption of this theory for the growth of crystals in the system occur according to the adjacent re-entry model. However, not all crystals in the system show this model of chains forming regular folds, as some crystals may also be connected by tie-chains and loops

• The initial atomic smooth face is not correct, as the surface of a real polymer crystal is irregular Over the years other models have been proposed to describe the polymer crystallization process, such as the multi-stage model by Strobl [35]. According to this model, the chain segments in the melt organize in a thin layer of mesophase or granular structure to the growth front. Later, this layer thickens until a critical size is reached. With an introduction of an intermediate stage for crystallization, the multi-stage model is very different from the HL-model. Another crystallization model was proposed by Muthukumar [34]. This model suggests that some chain segments contact the growth front, while the remaining parts of the chain quickly attach to the crystal surface. This means that there is a distinct diffusion for different parts of the chain. Moreover, there are other experiments and simulations which suggest the existence of some type of precursor before the formation of crystallites [31].

As mentioned earlier, crystallization can also be induced by stretching. When a coiled polymer is stretched, the number of accessible conformations reduces, which is often referred to as the origin of strain-induced crystallization (SIC) [36]. Flory [17] proposed a theory to describe SIC for a cross-linked network with high molecular weight, like vulcanized natural rubber. When the polymer network is stretched, some sections of the chain between the cross-link sites would also stretch, thereby reducing the conformation entropy. The assumptions of Flory's theory are:

• The system is in equilibrium • The crystallization does not happen before chains are fully stretched. Here, a chain is defined as a portion of a network extending from one cross-link site to the next one

• The crystallites align along the stretch direction

• Chains do not re-enter a crystallite by folding

• The conformations of amorphous chains follow Gaussian distribution

Crystal nuclei can be formed at a small strain, without much change in entropy. A crystal nucleus is defined as the temporary formation of a close-packed arrangement of segments or portion of chains in a plane perpendicular to the stretching direction. These nuclei can also act like giant cross-links (pivot points), that can occur at random locations in the network. When the system is further elongated, the conformation entropy will be less compared to the situation with no cross-links. To compute the conformation entropy change with respect to the fully crystalline polymer, Flory considers two hypothetical steps:

• Among n segments, the melting occurs for only nχ segments per chain and the chains ends are assumed to be free (i.e., not attached to any cross-links)

• The chain ends are assigned to the cross-link location within the elongated polymer

The entropy change for the first step is given by

S a = σ (n -χ) s f , (1.3)
where s f is the entropy of fusion per segment and σ is the total number of chains considered.

According to the lattice theory of thermodynamic properties of polymers

s f = k B [ln (γ -1) -1] , (1.4)
where k B is the Boltzmann constant and γ is the coordination number of the lattice. It represents the entropy of fusions arising from the randomness of the arranged segments. However, this does not include the contribution from internal changes that occur within the segments. The change in entropy due to the second step comes from the transformation of the amorphous chains to the elongated state. This theory predicts that the stretched rubber chains (reducing the possible chain conformations) results in an increased melting point, and thus favors crystallization.

Simulation studies

As mentioned earlier, polymer crystallization is a multi-scale process which can be studied using the computer simulation techniques like molecular dynamics (MD) and Monte-Carlo (MC).

The objective of MC simulation is to generate an ensemble of representative configurations under specific thermodynamic conditions for a given system. Consequently, the relevant thermodynamic observables, such as the free energy, can be calculated [37]. In MC simulations time is not defined and hence dynamical properties are difficult to access.

The experimental findings revealed the presence of macroscopically round and microscopically rough crystal surfaces for polyethylene [START_REF] Gedde | Fundamental Polymer Science[END_REF]. Based on these experimental observations, Sadler-Gilmer (SG) [38] put forward a crystallization theory by means of MC simulations of small molecules. A schematic representation a chain that attaches at three different positions which are labeled 1, 2 and 3 is shown in Fig. 1.11. Attachment 1 represents a stable crystal as its crystal thickness is greater than the minimum crystal thickness, whereas 2 and 3 attachments are not stable. The amorphous chain connecting to the attachment 2 are taut which prohibits the crystal thickness of this attachment. Attachment 3 can possibly grow on the lower side. The summary of SG-theory: when a chain attaches to a rough crystal surface, it does not choose suitable conformation for later stages of growth. The stem length may lie below the critical length (like attachments 2 and 3 in Fig. 1.11) where the chain folding and loops can prohibit further growth, like for attachment 2. These constraints forces the chain segments to find the suitable conformations for incorporation into the body of the crystal [START_REF] Gedde | Fundamental Polymer Science[END_REF].

The number of possible chain conformations leads to an entropic nucleation barrier, this is the essence of SG-theory [START_REF] Gedde | Fundamental Polymer Science[END_REF]. However, the general effects of neighboring lamellae and rearrangement of the chains in crystalline regions were not considered.

MD is a simulation method that solves Newton's equations of motion for interacting many-particle systems. The time evolution can be accessed via the trajectory of all particles in the system. A full-atomic simulation to study polymer crystallization requires extremely large computational resources and even then, it would be very difficult to simulate a simple realistic polymer. The main problem would be the disparity regarding timescales. Interesting phenomena such as a phase transition or phase separation in an experiment, happen on the timescale of a few minutes. However, a typical timestep in a MD simulation is in the range of 1 to 2 fs and so very long simulation times would be required to compare with experiments. Hence, coarse-grain (CG) methods have been used to access longer timescales. An example for a CG method is the united-atom (UA) model, where a carbon along with its attached hydrogen atoms are simulated as one bead and such UA simulation could have bigger timesteps compared to a full-atomic simulation. The idea of coarse-graining can be extended to mapping one or more monomers to a single bead, which typically would have even bigger timesteps compared to UA model. A schematic representation of CG models are shown in Fig. 1.12. Muthukumar et al. [START_REF] Muthukumar | Advances in Chemical Physics[END_REF] investigated crystallization of a polyethylene-like UA model, where a single chain exhibits chain-folding upon cooling. This represents the crystallization in the solution due to low density, as there is a lack of interaction between other polymers in the system. They reported that a single crystal based on a single polymer shows the lowest free energy in a folded state and free energy barriers dictates the initial lamellar thickness.

Muthukumar [START_REF] Muthukumar | Advances in Chemical Physics[END_REF] and Rutledge group [START_REF] Balijepalli | [END_REF] have investigate the initial stages of crystallization before the formation of stable crystalline domain. They observe two step process

-Objectives of this work

(two mesomorphic states) before forming stable crystals, Muthukumar labeled these as the baby nuclei and the smectic pearls.

Yi et al. [41] performed cooling of melts using linear polyethylene UA model for different chain lengths, to study the nucleation rate and the size of critical nuclei. They found that the nucleation rate is independent of the chain length but the growth rate decreases with an increase in the chain length. The initial thickness of the critical nucleus is reported to be shorter than the chain length where the crystal stems belong to different chains in the system. More details about the crystallization of polyethylene can be found in the review article [42].

Yamamoto [43], simulated UA polyethylene model to study fiber and network formation during stretch-induced crystallization. In the work specified by Yamamoto, the melt system was subjected to rapid elongation and then relaxed under isothermal conditions. Three stages of structural formation were observed, the first stage is critical nucleation, the second stage is the growth of lamellae, and the third stage is the formation of well-aligned stacked lamellae.

Hall et al. [44] used an existing Shinoda-DeVane-Klein (SDK) CG model [45] to study the structural properties of entangled and semi-crystalline polyethylene. Hall et al. demonstrated that the SDK model of polyethylene can be used to capture the polymer-scale structure of melt and crystalline phases. Their estimation for the melting temperature and thickness of lamella is in good agreement with the experimental findings. In this model phase transition occurs from a melt to a semi-crystalline state, but it does not exhibit the commonly observed crystal structure.

In recent years, new ideas and theoretical descriptions with respect to SIC have been put forward. Plagge and Klüppel introduced a relatively simple model describing the crystal formation and melting of stretched natural rubber, and in particular the dependence on temperature and cross-link density [16]. Sotta and Albouy took Flory's SIC theory and extend it successfully with a phase equilibrium type of description, in which they use the local draw ratio as the relevant order parameter [46]. Yokota and Kawakatsu studied nucleation and growth of nucleus in polymers with models based on classical nucleation theory and highlighted the significance of flexibility in polymers [47].

Objectives of this work

Natural rubber and its synthetic analog are widely used in tire manufacturing in both vulcanized and unvulcanized forms. For instance, when subjected to constant loading and unloading the tread bands of heavy-duty trucks harden the material and this hardening of rubber is believed to improve the mechanical properties [28]. In the case of thermoplastics like polyvinyl alcohol, industrial processing techniques like wet or melt spinning induces a phase transition which is also believed to increase the strength of the fibers [48]. This phenomenon of phase transition that occurs due to the application of some type of mechanical loading is referred to as strain-induced crystallization. The main objectives of this work are to study the SIC at the molecular scale by means of computer simulations using coarse-grained models and to characterize the induced structure.

Organization of this work

The organization of this work focuses on various aspects covered in our study. This work presents and covers the coarse-grained simulations for two different systems. As part of the literature study we have included the relevant experimental, theoretical and simulation works. We use bibliography as a reference database to validate our results and compare our measurements with the literature findings where-ever possible.

We start with the first Chapter by providing the general introduction of polymers and polymer crystallization. In the second Chapter 2 we have given an overview of two simulation techniques, molecular dynamics and dissipative particle dynamics. A brief description of different coarse-graining methods is also presented here.

In Chapter 3, we perform CG simulations of polyvinyl alcohol using a well established CG potential taken from literature. We investigate two modes of crystallization like temperatureinduced and strain-induced for different chain lengths. To quantify the induced structure we use various measurements like cluster analysis, local orientational order and so on. This is based on a recent article [49].

In Chapter 4, we perform atomic simulations and present a detailed method of developing a coarse-grained potential using the Statistical Trajectory Matching approach. The developed a CG model is used to perform continuous cooling and uni-axial straining at constant rates to investigate to two modes of crystallization. Similar to the previous chapter we also present various structural measurements.

In Chapter 5, we generate network configurations. Using the CG potentials obtained in the previous Chapter 4, we perform uni-axial deformation simulations using the dissipative particle dynamics technique.

In the last Chapter we summarize the entire work and present key findings and discuss some perspectives.

Chapter 2

Simulation Techniques

The topic of computer simulation techniques is vast, as there are variety of numerical methods that to model a system at different length and time scales. In the early days of computer simulations, a model of a liquid was simulated using the MC method by Metropolis in 1953 [50]. Over the years MC simulations have been used to model the systems that can be compared with the experimental data, like liquid argon.

A different simulation technique was introduced to obtain the dynamic properties of multi-particle system. One early attempt to study of dynamics of anharmonic, one dimensional crystal was performed by Fermi, Pasta and Ulam [51]. Later, in 1957 Alder and Wain-Wright [52] performed the first MD simulation to study the dynamics of a set of molecules using hard spheres, where the particles move at a constant velocity with perfect elastic collisions between each other. Several years later in late 1960's the properties of the Lennard-Jones model was investigated by Verlet [53]. The rise in efficiency of computational resources enabled the development of other simulation techniques, such as the finite element method (FEM) [START_REF] Taylor | The Finite Element Method for Solid and Structural Mechanics[END_REF], the elastic network model [START_REF] Munoz | [END_REF] and molecular simulations. Here, a choice has to be made based on the length and time scales of the relevant phenomena. In this work, we are interested in polymers at the molecular scale, therefore we focus on molecular methods like molecular dynamics and dissipative particle dynamics (DPD).

Molecular Dynamics

Classical MD simulation solves Newton's equations of motion for a group of interacting particles using effective forces. These forces or inter-atomic potentials can be determined by quantum mechanical calculations or using experimentally measured data. The equation of motion is expressed as

m i ¨ r i = f i , (2.1)
where m is the mass, r is the position and f is the force exerted on the particle i. The f i is calculated by differentiating the potential with respect to the position of a particle and is described as

f i = - ∂ U pot ∂ r i , (2.2)
where U pot is the potential energy function. The potential energy can include bonded and non-bonded interactions. The position and the velocities of all the particles in the system are computed at every step. There are several algorithms to solve Newton's equations of motion by updating positions and velocities like Verlet, velocity-Verlet and leap frog. All of them rely on the Taylor expansion of the particle positions at time t + δt, where δt is the short time integral referred as "time-step". In this work, we consider the velocity Verlet algorithm. The position of particles at time-step t + δt is computed using

r i (t + δt) = r i (t) + v i (t)δt + 1 2 a i (t)δt 2 , (2.3)
where v is the velocity and a is the acceleration of the particle i. The updated velocities are computed using

v i (t + δt) = v i (t) + δt 2 [ a i (t) + a i (t + δt)] . (2.4)
It should be noted that the new velocity v i (t + δt) can be computed only after computing the new position r i (t + δt). More details on these algorithms can be obtained in [START_REF] Frenkel | Understanding Molecular Simulation[END_REF].

Ensembles

MD simulations are used to study the natural time evolution of a classical system with N particles in a volume V. The quantities are time averaged during the simulation, which corresponds to the ensemble average in an ergodic system. Newton's equations of motion conserve the total energy E of a closed system, where N and V remain constant, such an ensemble is called a microcannonical or NVE ensemble. The original MD simulations were performed only in the NVE ensemble. However, it is desirable to simulate the system as close as possible to the experimental condition using ensembles like canonical (constant volume and temperature NVT) or isothermal-isobaric (constant pressure and temperature NPT) [START_REF] Binder | [END_REF]. The instantaneous temperature is related to the kinetic energy via the equipartition theorem

E k = N i=1 1 2 m i v 2 i = 3 2 Nk B T , (2.5)
where N is the total number of particles in the system, m and v are the mass and the velocity of particle i, respectively. In the ensembles where temperature control is required, a thermostat is introduced to maintain the average temperature of a system. There are many algorithms to control the temperature such as Andersen, Langevin (stochastic), Berendsen and Nosé-Hoover thermostats. The idea is to couple the system to a large heat bath.

In the case of an Andersen thermostat, heat bath and the system are coupled by allowing random collisions with imaginary particles. This is done by assigning a new velocity to a random particle drawn from a Boltzmann distribution with the desired temperature [58].

In the Langevin thermostat, two terms are added to the equation of motion mentioned in Eq. 2.1, a friction term that damps on the dynamics, and a noise term that allows to siphon the thermal energy into the system. The modified equation of motion is formulated as

m i ¨ r i = f i -m i γ i v i + W i , (2.6)
where γ i is the friction coefficient and W i is a stochastic fluctuating force through which the thermal energy transfer takes place. The thermal energy and friction are linked through the fluctuation-dissipation theorem to maintain the system at a desired temperature (see [37] for detailed description). Another way to control the temperature is to re-scale the velocity like in the Berendsen thermostat [59]. The velocities are scaled such that the rate of change of the temperature is proportional to the difference between the desired temperature and the actual temperature, expressed as

dT d t = 1 τ (T 0 -T) , (2.7)
where τ is a time constant that determines how the heat-bath and the system are coupled, T is the actual temperature of the system and T 0 is the imposed temperature. To re-scale the velocities, Berendsen proposed a modified Langevin equation of motion without random force

m i ri = f i + m i v i 2τ T 0 T -1 . (2.8)
Thus, the change in temperature is exactly equal to (T 0 -T) δt τ at every step and the velocity scaling factor is given by

λ = 1 + δt τ T 0 T -1 , (2.9)
It is important to mention that in the limit τ → δt the velocities are scaled almost instantaneously as the thermal energy exchange between the heat-bath and the system is large, whereas when τ → ∞ the thermal energy exchange is minimal, which means it would take very long time to achieve the desired temperature. Therefore, the total energy of the system remains nearly unchanged resembling the NVE ensemble [60].

Another method to perform constant temperature simulation is using the Nosé-Hoover thermostat. An additional degree of freedom s, which is a dynamic variable, is introduced to extend the real system. When the dynamic variable is greater than zero, it plays the role of a time-scaling parameter and the time interval between the real and the extended systems are expressed as d t = s -1 ( t) d t [61]. The thermal energy exchange between the extended system and the heat-bath is expressed via velocity scaling of the particles [62]. The relations between the real and the extended systems are given as

r = r , (2.10a) ṽ = ˙ r s , (2.10b)
where ṽ and ˙ r are the velocities of the extended and real systems, respectively. The coordinates r are the same for both systems. A modified Lagrangian for the extended system as postulated by Nosé is

L = N i m i 2 s2 ṽ2 i -U( r) + Q 2 ṡ2 -gk B T 0 ln s , (2.11)
where the first and second terms are associated to the kinetic and potential energies of the system, respectively. Note that for the first two terms velocities are multiplied by s, so they represent the energies of the real system. The third and fourth terms are the kinetic energy and the potential energy associated to s, respectively. The parameter Q determines the time scale of the temperature fluctuations and has the units of energy • time 2 . The factor g indicates the number of degrees of freedom, k B is the Boltzmann constant and T 0 is the imposed temperature. The equations of motion are derived from the Lagrangian Eq. 2.11

[63]

¨ ri = - 1 m i s2 ∂ U ∂ ri - 2 s ṡ ṽi , (2.12)
and

s = Q -1 s -1 N i m i s2 ṽ2 i -gk B T 0 . (2.13)
It is important to mention that the above two second-order differential equations are discretized based on the time-step in the extended system. Hence, the trajectory of the real system has to be interpolated at evenly spaced points. To solve the time scaling issue between the extended and real systems, Nosé-Hoover formulated a Hamiltonian expression corresponding to the Lagrangian of Eq. 2.11 and reformulated the equation of motion in terms of variables associated to the real system [61]. The Lagrangian equations of motion Eqs. 2.12 and 2.13 were written in terms of γ, where γ = γ = ṡ is a friction coefficient. The friction parameter becomes a dynamical variable similar to s and is a time-derivative term, it is postulated that the thermostat will provide a slow exchange of kinetic energy compared to the Berendsen algorithm. However, if the Q parameter is large, then the canonical ensemble will be obtained after very long simulation time, whereas if the same parameter is very small, unwanted oscillations may occur [61].

To perform a constant pressure simulation, the NPT ensemble is used. Since the pressure remains constant, the volume fluctuates. The pressure P is defined as [59]

P = 2 3V (E k -Ξ) , (2.14)
where V is the volume of the simulation box, E k is the kinetic energy and Ξ is the internal virial for pair-additive potentials

Ξ = - 1 2 N i< j r i j • f i j , (2.15)
where r i j = r i -r j and f i j is the force on particle i due to particle j. Similar to the Berendsen thermostat, the pressure is controlled by coupling to a pressure-bath. An extra term is added to the equation of motion and the pressure change is expressed by

d P d t = 1 τ p (P 0 -P) , (2.16)
where P 0 is the desired pressure and τ p is a time constant that couples the system to the external piston or the pressure-bath. To account for the changing volume, a coordinate scaling factor is added and the change in volume is expressed by

V = 3αV , (2.17)
where α is a scaling factor and 3 is the number of dimensions. The change in pressure is related to the isothermal compressibility β and is given as

d P d t = - 1 β V dV d t = - 3α β . (2.18)
An expression for the scaling factor can be obtained using Eqs. 2.16 and 2.18. The scaling factor α represents a proportionality of both the box size and the coordinates of particles per time-step. It is worth noting that Berendsen applied this approach to both isotropic and anisotropic systems. The later case yields a scaling tensor. Other barostats also follow a similar concept of implementation as their temperature counterparts. For example:

The Nosé-Hoover barostat is implemented using the same principles of the Nosé-Hoover thermostat, which is by adding an extra degree of freedom [62].

Force fields

A force field is a collection of different interactions between particles. In the early days of molecular simulations, the interaction between atoms were modeled as hard spheres, where the particles move at constant velocity with elastic collisions. Later, more complicated interactions were used to simulate liquids, small rigid molecules and polymers with a strong repulsion at short distances and a long range attraction, like the van der Waals interaction [START_REF] Allen | Computer Simulation of Liquids[END_REF]. Such a potential can be described by a Lennard-Jones (LJ) interaction

U LJ (r) = 4ε σ r 12 - σ r 6 , (2.19)
where r is the distance between two particles, ε is the depth of the potential well also referred as the interaction strength, and σ is the particle distance at which the potential U LJ (r) is zero (generally called the van der Waals diameter). The minimum of the LJ potential is located at the distance r min = 2 1/6 σ. In simulations, the LJ potential can be truncated at a distance r ≥ r min and the energy is shifted upwards to zero, such a model is referred as truncated-LJ potential. The truncated LJ potential is computationally less expensive compared to the full LJ potential, but still covers important physical features of the system of interest [START_REF] Stephan | [END_REF]. When the LJ potential only includes repulsive forces with the distance truncated at r min , it is called WCA potential [66]. Other variants of the LJ potential are also used in simulations with both repulsive and attractive parts similar to the full LJ potential. For example, σ r 12 in the equation 2.19 can be replaced by σ r 9 when softer interactions are desired. The LJ potential depends on the distance between a pair of particles. The forces should be obtained considering all pairs of particles. This would be computationally very expensive and hence a cutoff distance is introduced, where only the pair particles within a cutoff distance is considered in the molecular simulations.

For systems with charged particles, the LJ potentials are not sufficient to represent the long-range interactions. A simple approach is to supplement the pair-potentials with the Coulomb potential .20) where the distance between a pair of particles r i j , ε is the dielectric constant, q i and q j are the charges of particles i and j, respectively. The Coulomb interactions are long ranged and hence solvers like Particle-Particle-Particle-Mesh are used. The long ranged interactions are approximated by mapping all the particles in the system to a grid and then a potential is obtained by doing the Fourier transform of the charge distribution. A detailed presentation of long ranged solvers are presented in [START_REF] Mcdoniel | International Supercomputing Conference[END_REF]. The intramolecular potentials include bond, angle and torsion interactions. The simplest variety of bond potential is the harmonic type

U coul (r) = q i q j 4πεr i j , ( 2
ɸ h i j k θ r ij
U bond (r) = 1 2 k bond (r -r 0 ) 2 , ( 2.21)
where k bond is a stiffness parameter, r is the distance between the bonded particles and r 0 is the equilibrium bond distance. The spring constant or stiffness parameter determines the rigidity of the bonds. Generally in our systems, the bond lengths do not fluctuate much around the equilibrium bond distance. If the simulation time-step is big, bond lengths can increase leading to very high bond energies and large bond vibrations. There are several other types of bond interactions which are used in polymer simulations like,

• Finite extensible nonlinear elastic (FENE) [68] • Rigid bonds, follow a constrained algorithm like SHAKE [69] or RATTLE [70] • Quartic bond potentials

The angle interactions are found for particles that are part of two consecutive bonds, shown as angle θ between three bonded particles j, i and h in Fig. 2.1. Such interactions can also be modeled using the harmonic potential as given in equation 2.21, but (rr 0 ) 2 is replaced with (θθ 0 ) 2 , where θ is the angle formed between three particles and θ 0 is the equilibrium angle. Similar to bond potentials, there are many types of angular interactions, like the cosine potential

U ang l e (θ ) = 1 2 k angle (cos θ -cos θ 0 ) 2 . (2.22)
In order to appropriately model the interactions, more complicated angle potentials can also be incorporated using a tabulated form.

To model the interactions of four particles, two types of torsion interactions can be considered like proper dihedral and improper dihedral. If the particles are connected by three consecutive bonds it is proper dihedral, else it is improper dihedral. The proper torsion or dihedral angle is defined as an angle between two planes, depicted as φ formed in a quadruple set of particles k, j, i and h in Fig. 2.1. Commonly preferred dihedral conformations are 0 • , ±60 • and 180 • which are labeled as cis, gauche and trans states. Each of these conformations will result in a different polymer behavior and the Fourier torsion potential is one among many used to model these interactions in polymer simulations

U t orsion (φ) = m i=1 k t orsion i [1 + cos (n i φ -d i )] , (2.23)
where k t orsion i and d i are constants which represent the strength of interaction and the phase angle at equilibrium, respectively. The term n i is a positive integer indicating the periodicity and φ is the torsion angle formed between two planes. The goal of most simulations is to predict various properties of the system of interest. The interaction potentials are the most important input, as the measured quantities of the simulated system are only as good as the interactions used to model them.

Dissipative Particle Dynamics

DPD simulations were introduced by Hoogerbrugge and Koelman in the early 1990's to study the hydrodynamic behavior of liquids that was not accessible using classical MD simulations [71,72]. The atomic systems have more degrees of freedom compared to CG system, which are not relevant to study the flow of liquids. Hence, a coarse grained approach was applied, where a DPD particle represents a group of atoms and the equations of motion for the DPD particles was governed by certain collision rules. Later, Español and Warren modified the original DPD model with a continuous function in which the particles interact with three forces, a conservative force, a dissipative force and a random force [73]. A simple pictorial illustration is shown in 2.2, where these forces are shown acting between two DPD particles. A conservative force is derived from a potential, the dissipative force uses a friction parameter to reduce the relative velocity between the two particles and a random force acts equally on both particles. These two forces together is referred as "pair-wise Brownian dashpot" [74].

The extension of the DPD technique to polymeric systems was done by Schlijper, using a bead-spring model for a polymer-solution system, where a bead is a DPD particle that represents a group of atoms [75,76]. The effect of the time-step on different integration schemes to solve Newton's equations of motion and the phase separation was studied by Groot and Warren [77]. Following this seminal work, several studies were published to study the phase separation in linear block copolymers, where the formation of lamellar, hexagonal rod and micelles have been observed [74,78]. Generally, the length and the timestep associated to the DPD simulations are greater than the MD simulations. Hence, DPD simulations are classified as a mesoscopic numerical technique. Although, the value Figure 2.2: Schematic representation of interactions between two DPD particles. The three forces depicted from top to bottom correspond to, a conservative force, a dissipative force as a dampener and a random force shown as some noise, respectively. The Brownian dashpot conserves momentum and the forces go to zero beyond a selected cutoff. Image modified from Ref. [74].

of the time-step will depend on the polymer system and the relevant phenomena under consideration, as a rule of thumb it is bigger than those in the MD simulations.

Similar to the MD method, the time evolution of a set of particles is governed by Newton's equations of motion

d r i d t = v i , m i d v i d t = f i , (2.24)
where m, r, v and f are the mass, the position, the velocity and the force of particle i. The force contains three terms as mentioned above and each term is additive

f i = j =i F i j = j =i F C i j + F D i j + F R i j , (2.25)
where the sum is over all the particles within a specified cutoff radius and F i j is the force acting on particle i due to the interaction with particle j. The conservative force F c in the original DPD setup is often expressed as

F c i j = a i j w c (r)r i j , (2.26)
where a i j is the maximum repulsion between particles i and j. The unit vector ri j = r i j /r, where r i j = r i -r j . The conservative weight function is formulated as

w c (r) = 1 -r i j /r c r < r c 0 r > r c . (2.27)
In the original DPD model a fluid system is simulated using a set of particles where each particle represents a group of small molecules [79]. Later, a modified version of DPD was developed to simulate vapor-liquid systems which is referred as many-body DPD or MDPD. The conservative force in MDPD contains both attractive and repulsive forces between two particles [79]. A similar idea is applied in DPD simulation of polymers, where a polymer can be represented by a chain of particles connected by springs. The interaction between particles or beads is modeled using some type of bonded potential and the repulsive conservative force is switched off.

F D i j in Eq. 2.25 is a dissipative force, it depends on the relative positions and the velocities of particles i and j

F D i j = -γw D (r) ri j • v i j ri j , (2.28)
where γ is a friction coefficient, v i j = v i -v j is the relative velocity of particle i to particle j and w D (r) is a weight function depending on the distance r. The F R i j in Eq. 2.25 is a random force, which represents the effects of thermal fluctuations and depends on the relative position of the particles

F R i j = σw R (r)θ i j ri j , (2.29)
where σ = 2k B Tγ is a coefficient that represents the strength of the random force, k B is the Boltzmann constant and T is the specified temperature. The θ i j is a random variable with Gaussian distribution and unit variance, and w R (r) is the distance dependent weight function. Español and Warren [73] proposed a relation between the weight functions of the dissipative and the random force, and also used a simple model to determine the friction coefficient. However, other methods have been subsequently developed to obtain the friction coefficient and the CG potentials like a Bayesian method [80]. Groot and Warren [77] simplified the Español and Warren [73] method by using the same weight function for both random force and the conservative force

w D (r) = w R (r) 2 = (1 -r/r c ) s r < r c 0 r ≥ r c , (2.30)
where r c is the cut-off distance of the dissipative and random forces. It is common practice to use the same cut-off distance for all forces and s = 2. It was reported that the different value of s can modify the dynamical properties of a DPD system [79].

If forces between particles depend on the relative velocities, the standard velocity Verlet algorithm can not be used directly. Hence, Groot and Warren proposed a modified version of velocity Verlet algorithm [77]

r i (t + δt) = r i (t) + v i (t) δt + 1 2m i f i (t) δt 2 , (2.31a) ˜ v i (t + δt) = v i (t) + λ m i f i (t) δt , (2.31b
)

f i (t + δt) = f i r i (t + δt) , ˜ v i (t + δt) , (2.31c
)

v i (t + δt) = v i (t) + δt 2m i f i (t) + f i (t + δt) . (2.31d)
If the force is independent of velocity, the standard velocity Verlet algorithm can be recovered using λ = 1/2, where λ accounts for the effects of stochastic interactions. The predicted velocity of the particle i at time t + δt is denoted by ˜ v i (t + δt). In this modified velocity Verlet algorithm, the predicted velocity is used to obtain the force and then the velocity is corrected at the last step for each particle.

In this study we obtain the non-bonded conservative force, the parameters for bonded interactions and the friction coefficient using the Bayesian approximation referred as statistical trajectory matching [80] (more details will be presented in the following section). Note the different ensembles mentioned in sub-section 2.1.1 discussed above applies to the DPD simulations as well. However, exclusive thermostat algorithms are not required in the case of DPD technique as the dissipative force and the random force together act as thermostat.

Coarse-Graining Methods

The structural and dynamical properties of polymeric materials occur at different spatial and temporal scales. Hence, it becomes necessary to adopt multi-scale simulation techniques to model these materials. Consequently, a rigorous link has to be established between the chemical constituents composed of atoms and the relevant macroscopic mechanical properties [81]. However, setting up such a rigorous connection is still an ongoing work. The main difficulty originates due to the multi-scale nature of the polymer systems. The typical vibrations of covalent bonds happen at the time scale of less than picoseconds and at the length scale of an Ångström [81]. The relevant dynamics of a monomer typically occur in the range of tens of picoseconds and the size of a monomer is generally considered to be about a few Ångströms [81]. The size of a typical polymer network with entangled chains is reported to be on the length scale of millimeters to centimeters and the relaxation, and aging of these bulk systems occur in the range of seconds to even years [81]. Since the bulk behavior of polymeric systems depends on the individual polymers which in turn depends on its basic constituents (atoms), a multi-scale modelling approach becomes inevitable. Here we focus on the molecular modelling of polymers. Over the past couple of decades, various computational techniques like molecular dynamics (Section 2.1), dissipative particle dynamics (Section 2.2) and Monte-Carlo have been developed to model polymer systems using CG approaches.

The development of CG models relies on either a top-down or a bottom-up approach. In the top-down approach, the derived parameters are used to reproduce the experimentally measured properties of the system. A popular method in this kind of approach is the MARTINI force field [82,83]. This model was originally developed for lipids [83] and later extended to other types of systems [START_REF] Drobny | Handbook of thermoplastic elastomers[END_REF][START_REF] Marrink | [END_REF]. The bond and angle interactions are given by weak harmonic potentials and the non-bonded interactions are modeled using a shifted LJ potential. In addition, if charges are present in the system, a shifted Coulombic potential is used to mimic the effect of distance-dependent screening. In a MARTINI model, the CG beads or CG sites correspond to three to four heavy atoms with the exception of ring-like molecules. Although this approach might be able to reproduce and predict the experimental results, it has its limitations. A top-down model is limited by the available experimental data, because this model cannot reproduce the properties more accurately than the available information. In addition, different models may be capable of predicting the reference experimentally measured properties of the system, which therefore requires to measure other quantities that were not used to model the interactions and compare with the reference measurements to select the better model.

The alternative approach to obtain CG potentials is the bottom-up methodology where these potentials are derived from atomic or high resolution simulations. There are several bottom-up approaches, like the simplex parametrization [86], the Iterative Boltzmann Method (IBM) [87,88], the Inverse Monte Carlo method (IMC) [START_REF]Coarse-Grained Modeling of Biomolecules[END_REF], the force matching approach [START_REF] Ercolessi | [END_REF]91], the Blob model [92] and Statistical Trajectory Matching (STM) [80]. Apart from the above mentioned approaches, there are several other multi-scale computational techniques which are discussed in detail in the review article [81] and in books [START_REF]Coarse-Grained Modeling of Biomolecules[END_REF][START_REF]Coarse-Graining of Condensed Phase and Biomolecular Systems[END_REF]. A coarse-graining can be described in terms of mapping functions, that associates a group of atoms to a coarse-grained bead. The task of any CG method is to obtain a potential or a set of potentials that provide the same properties for both the atomic and the CG systems. In general, these set of target properties are obtained from a high resolution simulation. These properties can be for example, instantaneous forces, RDFs, structural distributions or pressure. The IBM, IMC and simplex parametrization are structure-based CG methods, where the target functions, like the RDF and other distributions, are replicated as good as possible. We can predict the viscoelastic properties from molecular level using the bottom-up approach We can design new polymeric materials for targeted viscoelasticity using the top-down approach Fig. 2.3 presents multi-scale hierarchical computational approaches developed to predict the viscoelastic properties of polymers. The force fields used in atomic simulations are generally obtained through ab initio methods like Density Functional Theory (DFT) and act as a bridge between the length scales of pico-and nano-meters. The coarse-grained potentials obtained from atomic or united-atom simulations bridges the length scales between nano-and micro-meters. Usually, the timestep considered to perform these simulations ranges from 1 fs to 50 fs, where the CG simulations are done with the timestep longer than 5 fs. The timestep mentioned above can only be treated as a thumb rule, it depends on factors like the system, phenomena of interest and the simulation procedure followed. A CG model can be used to compute the entanglement length by performing a Primitive Path (PP) analysis using a program like Z1 [95]. This PP network can further be mapped onto a tube model, which in turn bridges the length scale between micro-and milli-meters. The tube model can be used to derive a constitutive model for polymers at the macro-scale [81] with parameters obtained through CG simulations.

Simplex Method

The simplex algorithm [START_REF] Press | Numer. Recipes C[END_REF] is a multi-dimensional optimization scheme which was first used to optimize atomic simulation models for experimental data [86]. In the case of the automatic simplex algorithm a single-valued function is constructed f (p 1 ...p d ) to evaluate the set of force field parameters p 1 ...p d . The choice for f is the square root of the weighted sum of relative squared deviations [86]

f (p n ) = i w i 1 - P i (p n ) P i,t ar get 2 , (2.32)
where P i,t ar g et is the target value of the desired property P i that is obtained from experiments. w i is a weight to account for the fact that some properties, like density, may be easier to reproduce than others. The initial configuration should be as close as possible to the desired state, for example, the density of the system must be near the desired average value. If not, this algorithm might take a very long time to find the set of parameters that provide the appropriate values. The target function is first evaluated before starting the simplex algorithm. Then a simulation is performed using an initial guess of potentials. If the system is equilibrated according to some criteria, like the convergence of the density, then the target function value is evaluated again. It is an iterative process until the selected set of properties converges to the target values. In 2001, Faller et al. [START_REF] Faller | [END_REF] developed a coarse-grained model of trans-PI using the simplex method. To optimize the CG potential using the simplex algorithm, they define a single valued function f similar to Eq. 2.32, but use the RDF as a target function g t ar get (r) and which is obtained from an atomic simulation. The function to be minimized is the integral over the squared difference in RDF [98],

f = w(r) g(r) -g t ar get (r) 2 d r , (2.33)
where w(r) is a weight function which penalizes the deviations at small separation [99], g(r) is the radial distribution function measured at each iteration and g t ar get (r) is the target RDF obtained from a high resolution simulation. A drawback of the simplex technique is that a relatively small set of parameters, between 4 and 6, should be used to describe the optimized potential. A larger number of force field parameters increases the dimensionality for the simplex algorithm which increases the computational cost. Although, this method has been used to develop coarse-grained potentials of polyisoperene, much experience and physical intuition is needed to obtain the final set of potentials. In a comparative study performed by Reith et al. [87] it is reported that the quality of the potentials is not as good as that of IBM.

Iterative Boltzmann Method

The iterative Boltzmann method was developed to overcome the problems encountered with the simplex method [87,100]. However, similar to the simplex parameterization method, a structural measurement like the g(r) is used to obtain an optimized CG potential. The idea is to speed up the convergence of a desired function and at the same time get rid off the limited number of force field parameters as imposed in the simplex method [100]. In the IBM method, the initial guess of the potential is done by Boltzmann inversion of the target radial distribution function which is obtained from an atomic simulation. The potential obtained from inverting the target g(r) between the non-bonded monomers accounts for the local structure. It has been shown that simple packing arguments can be used to explain the local orientation correlations in dense melts [98]. However, the pairwise potential obtained from a simple Boltzmann inversion of the target RDF is incorrect for a bulk system, thus this is used as an initial guess.

To derive the CG potential using IBM, a polymer system is simulated using an atomic model. The atomic configuration is mapped onto a coarse-grained system and the probability density functions are measured for bond length, angle and RDF. The measured probability density functions are Boltzmann inverted to get the initial guess for the potential to start the CG simulation. Boltzmann inversion is performed using the following equation [100]

U(R N ) = -k B T ln p R (R N ) , (2.34)
where U(R N ) is the potential energy function, R N are the coordinates and p R (R N ) is the probability distribution of the N CG beads. In most cases the p R depends on four variable sets which are 3N dimensional functions of the coordinates of the CG beads. The four variable sets correspond to pair distance r, bond length l, angle θ and dihedral angle φ. It is assumed in IBM that these four variables are independent of each other and the potential function of the CG model becomes

U(R N ) = U(r, l, θ , φ) = -k B T ln p R (r, l, θ , φ) , (2.35a) U(r, l, θ , φ) = U(r) + U(l) + U(θ ) + U(φ) . (2.35b)
As the idea is to reproduce the distribution function of the atomic model as accurately as possible using the CG model, an iterative process is followed

U n+1 (q) = U n + ∆U n (q) , (2.36a) ∆U n (q) = k B T ln p n R (q) p t ar get R (q) , (2.36b)
where q = r, l, θ , φ and p t ar g et R

is the target distribution function computed from the atomic simulation. The distribution functions p n R need to converge to the target distribution functions, after some iterations. Here n is the iteration step.

A typical workflow followed in the IBM can be summarized as follows [81]. The target distributions are obtained from the atomic simulations. These distributions are Boltzmann inverted to get an initial guess for the potential functions of the CG model. Then canonical MD or MC simulations are performed to generate the trajectories of the CG model. From these trajectories the distribution functions p n R are calculated and the potentials are updated using the Eq. 2.36. If the convergence criteria is fulfilled the iteration process is stopped and the optimized CG potential is obtained. If not the algorithm proceeds to the next iteration step. Though this method is widely used, the potentials derived may not be unique and need additional fine tuning to obtain the desired pressure. It is important to note that the pressure obtained in the canonical ensemble simulation using the CG potentials is often different from that of the atomic simulation. This is generally the case as the pressure is not optimized during the IBM iteration. In order to solve the pressure discrepancy, a correction term is added to the non-bonded potential. This is also an iterative procedure which is done by a second set of simulations. This iteration is continued until the relative difference in densities between atomic and CG systems is smaller than a threshold value.

It should be taken into account that the additional pressure correction term included in the pair potential does not change the measured bonded distributions with respect to the reference distributions.

Inverse Monte Carlo Method

Another structure-based method is the inverse Monte Carlo or inverse Newton method. Similar to the IBM, the goal of this method is to obtain the CG potentials that can reproduce the structural measurements like bond length, angle and g(r) as close as possible compared to that of the atomic simulations. In the IMC the relation between the set of parameters and the averages is given by a non-linear multi-dimensional equation. To find the set of optimized parameters from the set of averages obtained from MD or MC simulations is referred as the inverse problem. The set of parameters are iteratively obtained using a Newton-Raphson method.

The CG potentials for polymer systems are optimized using the IMC method where the interactions are in the form of tabulated potentials. These tabulated potentials can be described as constant grids. The total potential energy of the system is given by

U total = α E α A α , (2.37)
where the potential at each bin, E α , corresponds to a bonded or non-bonded energy with the corresponding occurrence A α . A tabulated pair potential is generally expressed as a step-wise function of distance [START_REF]Coarse-Grained Modeling of Biomolecules[END_REF], but may also include an interpolation scheme. The initial estimate for the CG potentials are obtained using the Boltzmann inversion of the probability distributions of bond length, angle, dihedral and RDF from the reference trajectory. Using this initial set of force field parameters a CG simulation is performed. The deviation ∆ 〈A α 〉 of the average properties are computed at each iteration step n using

∆ 〈A α 〉 = A n α -A t ar get α . (2.38)
Next, the update ∆E γ is obtained after using the Newton-Raphson method to solve [89]

∆ 〈A α 〉 = γ ∂ 〈A α 〉 ∂ E γ ∆E γ + O(∆E 2 γ ) , (2.39)
where

∂ 〈A α 〉 ∂ E γ is the matrix of derivatives. It is assumed that γ = 1, ...,
M is finite number of potential parameters and α = 1, ..., M is the finite number of target properties. Note that A is a concatenated vector with populations sampled at every bin of the tabulated potential which includes bonded and pair terms. The vector ∆E contains iteration updates to all tabulated terms [101]. The Jacobian is expressed in terms of canonical averages which can be evaluated by any simulation method with a given set of force field parameters. The Jacobian is calculated using

∂ 〈A α 〉 ∂ E γ = - 〈A α A γ 〉 -〈A α 〉〈A γ 〉 k B T , ( 2.40)
where the indices α, γ represents all the interaction types (bonded and non-bonded). These interaction types include all CG beads (types of beads) and for each set of CG types, covers the relevant range of distances or angles [START_REF]Coarse-Grained Modeling of Biomolecules[END_REF]. The average values of 〈A α 〉 and crosscorrelation 〈A α A γ 〉 can be obtained from a CG simulation. The obtained ∆E is used to update the potentials. This procedure is followed until a convergence criteria is satisfied.

Other CG methods

The choice of a CG mapping is not unique. It is possible to define a CG bead that represents several monomers of a polymer. Such a model with a higher level of coarse-graining can extend both the time and the length scales of the CG simulation. The blob model is one such systematic CG approach. In this model the CG bead is considered to be a spherical blob, which represents the center of mass of several consecutive monomers of a polymer. All the potentials for this model are derived systematically from atomic simulations using the IBM and the IMC method.

In general, the motion of an entangled polymer is confined by adjacent chains [102]. In the literature, this restricted environment of a polymer is referred as a "tube". If an arbitrarily large number of monomers is considered to be a blob, then the size of the blob can exceed the diameter of the tube, which therefore should not be mapped as one CG bead. Since the idea of the blob model is to simulate the long-time relaxation behavior of a polymer system, which is achieved by considering large integration timestep, the size of the CG-bead should be large enough to allow the bigger timestep [81].

Padding and Briels [92, 103] modeled a polyethylene system with each CG bead representing about one-third of the entanglement length (each bead has 20 monomers). Since a large number of monomers is mapped as one CG bead, the dihedral interactions between these CG beads are weak and thus ignored in the simulation. The bond interactions include both a repulsive and an attractive term and the bonds become very soft due to high level of coarse-graining and thus unphysical bond crossings may occur. In order to overcome this problem the TWENTANGLEMENT algorithm was developed. More details of this algorithm are available in [92].

The earlier studies using the DPD technique did not apply any constraint to prevent chain crossing, therefore DPD could not be used to study entangled systems [81]. Nikunen et al. [104] developed a simple and computationally efficient method to avoid chain crossing in DPD simulations. Goujon et al. [105] developed a repulsive potential between segments to avoid bond crossing, and later [106] employed anisotropic potentials in DPD simulations for a linear polymer system. In that study, the shapes of the beads is chosen to be ellipsoidal based on the parallel and perpendicular components of the radius of gyration tensor.

However, the volume of a bead is fixed. The anisotropy factor h, which is defined as the ratio between the parallel to the perpendicular components of the gyration tensor, forms an additional input parameter for the simulation.

Statistical Trajectory Matching

Most of the bottom-up approaches discussed above, like the IBM and the IMC, use an iterative fitting procedure to obtain the optimized CG potentials based on the structural measurements of the high resolution simulation. These CG potentials are used to perform either MD or MC simulations. Dequidt and Canchaya [80] introduced a new bottom-up method for optimizing the CG potentials for DPD simulations which is referred to as statistical trajectory matching. The idea of STM is to match the whole reference trajectory, which is obtained from atomic or UA simulations, using the Bayesian method by Español and Zùñiga [107]. The idea of using Bayesian inference to optimize the CG model was first suggested by Liu et al. [108]. The reference atomic trajectory is mapped onto a coarse-grained trajectory that contains the coordinates of the CG beads. STM aims at statistically reproducing the total force on each CG bead during a timestep.

A DPD potential has three force terms, the conservative, the dissipative and the random forces as discussed in Chapter 2. It is important to point out that the CG potentials obtained from STM do not have the same form as that of the original DPD model introduced by Warren-Groot [77]. This will become evident from the description below. The basic principle of STM is Bayes' theorem. It states that the probability to find the set of parameters x, knowing the reference trajectory T r e f = {R 0 , R 1 , ....} is proportional to the probability of this trajectory T r e f being reproduced by a model using these set of parameters x, given the initial conditions R 0 and R 1 [80],

P(x|T r e f ) ∝ P(R 2 , R 3 , ...|x, R 0 , R 1 ) , (2.41)
where R is the set of coordinates of the CG beads and the indices 0, 1, ... represent the consecutive configurations in the reference trajectory. The time between the successive configurations in the reference trajectory will be the timestep that should be used in CG simulation after obtaining the optimized potentials. Note that the position R and velocity V vectors are 3N dimensional. The optimal set of potential parameters is the one that maximizes the likelihood of P(x|T r e f ).

Using the definition of conditional probabilities [80],

P(R 2 , R 3 , ...|x, R 0 , R 1 ) = P(R 2 |x, R 0 , R 1 ) × P(R 3 , ...|x, R 0 , R 1 , R 2 ) . (2.42)
The velocities and forces are computed using a first order time integration algorithm like velocity-Verlet or Leap-Frog. It is worth mentioning that technically any other time integration schemes like the modified velocity-Verlet method introduced by [77] as discussed in Chapter 2 can also be used. However, a systematic study about the dependence of STM potentials on different time integration scheme is yet to be done. The force or acceleration acting on the particles at time t can be deduced using the position of particles in successive configurations of the reference trajectory as expressed by,

F t = M R t+1 -2R t + R t-1 δt 2 , ( 2.43)
where M is a diagonal matrix that contains masses of the CG beads and δt is the difference in the time between successive configuration, i.e., it corresponds to the time between the configurations saved during the MD simulation. The position of the particles R 2 can be substituted by the force F 1 using Eq. 2.43 in the probability function mentioned in Eq. 2.42.

Then the probability of finding the position R 2 , when the R 0 and R 1 are known, can be expressed as [80],

P(R 2 |x, R 0 , R 1 ) ∝ P(F 1 |x, R 0 , R 1 ) = P(F 1 |x, R 1 , V 1 ) , (2.44)
where V 1 is the velocity of the particles at a particular discrete time. This is an affine transformation where the Jacobian is a constant. Similar to any DPD model, it is assumed that the force F t depends only on the positions R t , velocities V t and random numbers with no time correlations [80]

P(x|T r e f ) ∝ t≥1 P(F t |x, R t , V t ) .
(2.45)

If the positions, velocities and the force field parameters are known, the force between the CG beads are also known. However, the random forces are unknown and the probability function to be maximized is given by [80]

P(x) ∝ t P(F R t = F t -F C t -F D t |x, R t , V t ) , (2.46)
where F C t , F D t and F R t are the conservative, dissipative and random forces, respectively. In STM it is assumed that F D t is linear in the velocity [80]

F D t = -Γ t V t , (2.47)
where Γ t is a symmetric friction matrix that depends on the positions of the coarse-grain particles. The fluctuation-dissipation theorem requires that the random forces satisfy the following condition

[73] 〈F R t F R * t 〉 = 2k B T Γ t δ(t -t ) , (2.48)
where * indicates the transposed vector. The random forces are assumed to follow a normal distribution [80]. The probability that a model reproduces the trajectory exactly can be expressed as [80]

P(x) ∝ t 1 det(2πΓ t ) exp - δt 4k B T F R * t Γ -1 t F R t .
(2.49)

Maximizing P(x) amounts to minimizing L(x) = -log P(x). Note that the volume can be included in the set of CG parameters if the reference trajectory is obtained in the NPT ensemble. This is more likely to reproduce the exact reference density. The analytical solution for the minimizing of L(x) is presented in [80]. The function L(x) is given as

L(x) = 1 2 t log detΓ t + δtF R * t Γ -1 t F R t 2k B T , (2.50)
where δt is the timestep used in the CG simulation,

The STM approach can also be extended to obtain temperature dependent CG potentials. In such a case the reference trajectories are obtained at different temperatures and subjected to the same fitting procedure as described above. More information about this can be found in Ref. [109].

The workflow of STM is relatively straightforward. First, the MD simulations are performed using an all-atom or UA model, and MD configurations are saved at a regular interval. The time-interval at which the MD configurations are stored in the MD trajectory will be the timestep used for the CG simulations. Next, a CG mapping scheme is applied to the MD trajectory. This CG mapped trajectory, which is also referred as the reference trajectory, is subjected to the STM fitting process to obtain the CG potentials. The optimized CG force field is then used to perform DPD simulations.

The STM approach has been used to study different systems. In the study performed by Solano Canchaya et al. [109] a simple molecule, n-pentane, is considered to calculate the bulk and interface properties. The goal of this study was to demonstrate the transferability of CG potentials across thermodynamic conditions. For this purpose, they developed temperature dependent and temperature independent models, where they found that the temperature dependent model performs better at reproducing various properties when compared to reference MD model. Kempfer et al. [110] applied STM to model cis-polybutadiene at different temperatures and pressures. The thermal expansion coefficient and isothermal compressibility were well reproduced with respect to the MD and the experimental values. STM was also used to study the impact of different levels of coarse-graining, where the level of coarse-graining ranged from one monomer per bead to 5 monomers per bead. It is reported that for the low level coarse-graining the thermodynamic, structural and dynamical properties are well reproduced when compared to those of the reference MD computations [111]. The STM procedure is employed to obtain a CG model to investigate the structural properties of an interfacial region between silica nanoparticles and cis-PB [112]. The dynamics of the silica nanoparticle is well reproduced at short time scales in the case of low degrees of coarse-graining, i.e., one monomer per CG bead. The CG model obtained for the cis-polybutadiene system in Ref. [111] is used to investigate an interface model and examine the adsorption behavior [113]. For this purpose, the cis-PB in this interface model is confined between solid cuprous oxide layers.

Chapter 3

Crystallization of polymers by coarse-grained molecular dynamics simulations 3.1 Introduction

Polyvinyl alcohol is a hydrophilic, biodegradable, and atactic synthetic polymer that exhibits crystallinity. In the early 1920s PVA was discovered and since then it has been used in various industries, like food, textile, paper manufacturing, and in the production of polarizer films [114][115][116]. Mainly due to some important properties like nontoxicity, bioadhesive, significant uni-axial tensile, and surface characteristics [117].

It is reported that water is the main modifier for melt spinning of PVA and the drawing process enhances the mechanical properties in PVA fibers [48]. Wu et al. [48] suggest that during the drawing process, the chains orient along the fiber axis, and the structure changes with crystal orientation, lateral periodicity in the amorphous region, and crystallization, which accounts for the increased performance of processed fibers compared to as-spun fibers.

There are many experimental studies [48,[117][118][START_REF] Miyasaka | Structure in Polymers with Special Properties[END_REF][START_REF] Hong | [END_REF] concerning the structural evolution of PVA fibers and films. Miyazaki et al. [118] performed uni-axial stretching of PVA films with the thickness of 0.075 mm in water. They hypothesized that the PVA films consist of micro-fibrils, in which the crystalline and the amorphous regions stack in an alternating fashion, and the SIC occurs due to the extension of the amorphous regions, which results in increased stress with increasing strain.

Since the 2000s, molecular simulations have become a powerful tool for investigating molecular mechanisms occurring at the polymer chain scale in crystalline polymers [121]. However, some of these simulations use an atomistic resolution and therefore are limited by the time and length scales that can reasonably be sampled with the available computing infrastructures. The need to model long chains, large system-sizes and reasonable strain rates has sped up the development and use of coarse-grained models. With these lower resolution models, the basic element known as bead may represent one to several monomers depending on the degree of coarse-graining assigned to the CG model.

Furthermore, new industrial constraints related to recyclability, biodegradability, and bio-sourcing, raise many challenges to engineers, physicists, and chemists in their ability to synthesize or gauge new polymers exhibiting SIC with a relatively good ecological consequences. Finding new structures is vital from both sustainability and economic perspective.

In order to tackle such challenges, it is important to develop new materials at the microscopic level by retaining the most relevant aspect of the underlying chemical structures at the atomic level. This requires the ability to predict the crystallization phenomena of the new microscopic structures and its effect at the proper mesoscopic scale. It would be ideal to study any polymer system with all its fundamental makeup, but the constraint of time and length scales has influenced the development of different coarse-grained models. Although CG models have been used widely to access longer chain behavior, it has its own limitations. Olssen et al [122], performed MD simulations for tensile loading using atomic and UA-models of polyethylene. It was shown that during the tensile loading the yielding of semi-crystalline lamellar structures resulted in the formation of cavities in the case of atomic simulation, while the UA-model resulted in chain melting and transport. They concluded that, C H and H H interactions play a vital role and that the two models can result in different behaviors.

The CG model of PVA has been proposed in 2001 by Meyer et al. [123,124]. It is an excellent model system to investigate SIC due to its relative simple chemical nature and that can be considered on the atomistic level, in a united-atom approach, and on the coarsegrained level in the form of a bead-spring-model. Here we will focus on the coarse-grained model, because the processes we are interested in take place on the longer time-scales.

The model has been used to study the thermal induced crystallization of PVA and its dependencies on chain length [125]. On cooling the model exhibits the arisal of local structure that can be characterized in terms of folds and stems [124,125]. By exposing the polymers in simulations to a cyclic heating and cooling process at faster rates, the model can be used to explore the nucleation and growth of lamella [126,127]. Other topics of interest are the entanglement behavior [128,129], the extension to branched polymers [130], and various studies on the deformation of samples at constant temperature [131][132][133]. Here we aim to address fundamental questions concerning the response of PVA under a uni-axial strain at temperatures at and above the crystallization temperature. We compare results for short and long polymer chains in order to elucidate the impact of the polymer length.

The Model

We use a commonly used coarse-grain model that was developed for semi-flexible polyvinyl alcohol (PVA) by Meyer et. al [123,124]. The coarse-graining of PVA is based on an atomic simulation of 48 atactic decamers [CH 3 CHOH [CH 2 CHOH] 9 CH 3 ] in an orthorhombic simulation box with periodic boundary conditions using the YASP program [134]. The pressure and temperature is controlled by a Berendsen barostat and thermostat at 101.3kPa and 500K with a damping time of 5.0 and 0.2 ps, respectively, and the timestep used was 2 fs [99]. The atomic interactions used are modeled as rigid bonds using the SHAKE algorithm, a harmonic angle potential and a dihedral-angle potential. The non-bonded interactions are modeled by LJ potentials with a cut-off distance of 1.1nm [135].

In the CG system, one monomer corresponds to one CG-bead, where the mapping center is located on the backbone carbon of the CHOH group, as the resulting bond length distri-bution is independent of the angle distribution for this mapping center. The intra-chain distributions of the bond and angle of the atomic trajectories are used to find the suitable CG potentials by Boltzmann-inversion. A harmonic bond potential with the preferred bond length of 0.26nm and the stiffness parameter k bond = 2704k B T 0 /σ 2 was obtained, where k B is the Boltzmann constant and T 0 = 550K is the reference temperature. The tabulated angle potential is shown in Fig. 3.1 with local energy minima corresponding to the gauchegauche, the trans-gauche, and the global minimum for the trans-trans state at approximately 95°, 126°, and 180°angles, respectively. The tabulated angle potential effectively includes the torsion interactions found at the atomic level. The radial distribution function (RDF) of the atomic simulation is computed excluding the four nearest neighbors of all possible backbone carbons of the CHOH group and this RDF is treated as the target function. The CG potential is optimized to fit this target function, using a 6-9 Lennard-Jones potential, that is truncated and shifted to zero at the cutoff distance r c = 1.02σ

U LJ (r) = ε σ LJ r 9 - σ LJ r 6 , (3.1)
where σ LJ = 0.89σ. The interaction strength ε is fixed at 1.511k B T 0 . The chain diameter is chosen as the unit-length σ = 0.52nm, which is roughly the distance between two carboncarbon bonds in the backbone. The unit of time is chosen to be τ = σ m/k B T 0 = 1.61ps, and the temperature is given in reduced units T = T real /T 0 .

We perform all the simulations in this chapter using LAMMPS, a particle dynamics code [136] with a timestep of δt = 0.005τ and periodic boundary conditions. The temperature is controlled by a Nosé-Hoover thermostat with a damping-time of 100δt and the atmospheric pressure is imposed using a Berendsen barostat with a damping-time of 1000δt. Two systems of chain lengths of 100 and 1000 monomers are considered, where the short chain system has 2700 chains and the long chain system has 100 chains, respectively.

The initial configurations were generated as random polymers by taking into account the bond and angle interactions within each chain. These systems were equilibrated to remove possible artifacts. An example for an artifact is the octahedral structure formed by the beads of different chains which is schematically shown in Fig. 3.2. To eliminate such an artifact, soft repulsion potential is applied to begin the equilibration process and gradually increase the repulsion strength of the 6-9 potential. The short chain systems were generated with 100 polymers and then subjected to the above mentioned process. Later, these equilibrated configurations were periodically replicated to have 2700 and 21600 chains. These big systems were further subjected to longer equilibration to remove any existing correlations. The systems were equilibrated for 0.25 µs in the case of short chains and 25 µs for long chains at the temperature T = 1.

Temperature induced crystallization

Macroscopic scale

The equilibrated systems were subjected to different constant cooling rates Ṫ between 10 -2 τ -1 and 10 -6 τ -1 at atmospheric pressure. Note that even the slowest quench rate we have considered, is virtually instantaneous in real units. To understand the behavior of both the systems at different cooling rates, the measured volume per monomer as a function of the imposed temperature is shown in Fig. 3.3. To compare the results of short and long polymers, the average volume per monomer for the long polymer have been shifted up with a constant density 0.00269/σ 3 . The difference in density arises from the fact that the end monomers of a polymer form a single bond, whereas the internal monomers have two bonds. Hence, the short chain systems have a slightly higher volume per monomer in comparison with that of the long polymers. As a reference, the unmodified result of long polymers at the slowest cooling rate is also presented in the same graph. As expected for the fast cooling rates, the measured volume per monomer shows a linear dependency with respect to the imposed temperature and the systems become glassy. However, for the slower rates of cooling, the polymers have more time to relax and adjust with regard to their changing temperature thereby, exhibiting a phase transition to the semi-crystalline state, which is indicated by the drop in volume or equivalent increase in the density. For the short chain system, the results are averaged values over 5 independent realizations, whereas for the long polymers only a single set of quenches have been used. Only for the slowest cooling rate, the results for both the chain lengths do not match. The long polymers need more time to relax, while, the short polymers have enough mobility to relax. We estimate the crystal-lization temperature (T C ) from the location where the curves leave the linear domain. This T C depends on both the cooling rate as well as the polymer length, which are 0.81 and 0.82 (≈ 450K) for the short and long polymers, respectively. It is to be expected that for slower cooling rates the crystallization process will begin at slightly higher temperatures. To illustrate this prediction, we perform a cooling simulation using the quench rate of 10 -7 τ -1 for a long chain system and the resulting cooling curve is shown in Fig. 3.4. The observation that the long chain system starts to crystallize at higher temperatures is associated with their reduced mobility. It can be noticed that upon slow cooling the volume occupied per monomer reduces compared to that of the constant cooling rate of 10 -6 τ -1 , as the chains have more time to relax to form compact structures. Experimentally the melting temperature of PVA is estimated to be about 500K [137,138], which is underestimated by nearly 4% in our model as shown by heating curve in the same figure. We estimate the melting temperature to be about 480K from the location where the heating curve leaves the initial linear domain. Upon further heating at a constant rate, a sharp increase in the volume can be observed and the heating curve joins the cooling curve above the reduced temperature of 0.9.

Polymer scale

Snapshots of a long polymer during continuous slow cooling at different temperatures are shown in Fig. 3.5. The polymer has a coiled conformation in the melt with relatively few folds. Note, the formation of straight segments within a polymer in general (for this model) begins at T = 0.82 which is the transition temperature from melt to lamella. However, chain folding might not occur in every chain. It can be observed that the formation of straight segments for this particular polymer shown occurs at T = 0.80. Upon cooling the melt system at the slow-rate, the chain folding begins in some polymers which can become part small crystallites. In this particular polymer, the chain folding occurs below the T C . Upon further cooling, more straight segments occur, with different orientations, that are connected by disordered sections within the same polymer. To quantify the effective size of a polymer, the radius of gyration R g can be computed for each chain

R g = R 2 x x + R 2 y y + R 2 zz , (3.2)
where R x x , R y y and R zz are the diagonal elements of the gyration tensor. In Fig. 3.6 the evolution of R g for both chain systems as function of temperature is shown for different cooling rates. During the fast cooling rates, the polymers do not have sufficient time to relax and adjust to the changing temperatures. Hence, the polymers get frozen and form a glassy state. Upon slower cooling, the internal structures of the polymers are affected above their crystallization temperature, because they have more time to relax. To capture the structural changes within the polymers, we measure the bond length and angle distributions during cooling. In Fig. 3.7 the distributions of the slow cooling rate for a short chain system is shown. As the temperature reduces the bond distances tend to become narrower and the preference for the bigger angles increases. Hence, it appears that it is mainly the angle potential that is responsible for the formation of straight segments within a chain. The angle potential has a global minima at 180 • which results in a "transport" towards bigger angles, leading to more extended polymers.

It is worth noting that the average radius of gyration does not characterize the behavior of individual polymers due to the variation in sizes. In Fig. 3.8 the probability distribution of the radius of gyration during slow cooling rate for short chain systems is measured in temperature windows of ∆ = 0.04.

In the melt the R g distribution has a uni-modal profile with a range of values between 2σ and 9σ. On cooling the measured R g distribution still remains uni-modal, but the range of R g value changes compared to the melt. This uni-modal behavior can be seen for configurations well below the crystallization temperature. However, for the lowest temperature windows the shape of the distribution is not uni-modal and it is tempting to assume that this is merely due to poor statistics or configurations not being independent. The results presented here are obtained from 5 different realizations and the error-bars are shown in the same graph. We have also considered much bigger system of 21, 600 polymers and the results are the same with respect to 2, 700 chain systems. We have confirmed that the R g values of individual polymers remain dynamic down to temperatures T ≥ 0.75, but below this temperature the polymer conformation no longer de-correlates within the simulation time.

The shape of the probability distribution of the radius of gyration is associated to the polymer length that is cooled below the crystallization temperature. For the short polymers, the peaks close to 7σ and 5σ corresponds to well defined structures, i.e., single and double hair-pin, respectively. Examples of such conformations are shown as insets of Fig. 3.8. By analyzing the average R g we have implicitly assumed a spherical symmetry for the polymers, but this is in fact not true. Even in the melt the overall shapes which is characterized by the components of the gyration tensor, remains anisotropic and the aspect ratio increases on cooling. The relative anisotropy shape is defined as

κ 2 = 3 2   λ 4 1 + λ 4 2 + λ 4 3 λ 2 1 + λ 2 2 + λ 2 3 2   - 1 2 , ( 3.3)
where λ 2 1 , λ 2 2 and λ 2 3 are the eigenvalues of the gyration tensor and they are ordered as

λ 2 1 ≤ λ 2 2 ≤ λ 2 3 .
The value of relative anisotropy is bounded between zero and one. Zero corresponds to a spherical shape and one characterizes a straight segment like a rod. The measured relative anisotropy as a function of the imposed temperature for different cooling rates is shown in Fig. 3.9.

The change in average R g during cooling is less for long chain systems compared to the short polymers. Moreover, for the fast cooling rates shown in Fig. 3.6 it appears that the average R g value reduces with temperature for the chain length of 1000 monomers. For an Figure 3.9: The relative shape anisotropy for polymers of (a) 100 and (b) 1000 monomers during cooling. isolated polymer the bigger angles are preferred, which results in the polymer extension. However, in the case of a bulk system this is severely hindered due to the tube-confinement and the presence of entanglements. Luo and Sommer [125] have estimated the entanglement length to be between 20 and 30 monomers for a chain length of 1000 monomers. Our long chain system has only 100 polymers and we have considered only a single realization and thus it suffers from too poor statistics for an accurate representation of the average behavior. The long chain systems are also affected by significantly slower dynamics and hence would need much slower cooling rates, which are at present outside of our computational reach.

The preference for the bigger angles upon cooling, enhances the formation of short straight segments. These straight segments are called stems and are regarded as the building blocks of the lamellar structure [127]. The stems start to occur just above the crystallization temperature as shown in Fig. 3.10. To identify the stems, we use a simple criterium where a particle is labeled as being part of a stem whenever the surrounding segment exceeds a reference distance. For a given particle we consider two monomers that are N b beads away on either sides along the polymer. Since the average bond distance is b 0 , this gives a reference distance in terms of the maximum distance α (2N b b 0 ), where α is a proportionality coefficient. We have used the values of α = 0.9 and N b = 3 to compute the stem length. For the beads near the end of the chain, the segment is truncated at one side and the reference distance is modified accordingly. The consecutive beads along the chain that are identified as part of a stem belong to the same stem and a single polymer can have multiple stems.

We measure the number of stems and the stem length in each temperature window according to the above mentioned approach. In Fig. 3.10 the length-weighted average stem length is shown for both polymer lengths at different cooling rates. After reaching the crystallization temperature, the average stem length increases rapidly until the semi-crystalline regime is reached for the slow cooling rate. The average stem length for different cooling rates are similar for both the chain lengths. For slow cooling rates, the average stem length is approximately 20 units, which is in agreement with the expected value [127,131,139] and also comparable with the results for polyethylene [140]. According to Hoffman-Lauritzen nucleation theory [31,141], the stems longer than the minimum lamellar thickness will act as nucleation sites on which other stems are deposited and from which a lateral growth of the crystal is initiated.

Monomer scale

The radial distribution function or pair-correlation function g (r) of the monomers provides information about the local positional order. The results of the pair-correlation functions for a long chain system are computed at different temperatures during the slow rate of cooling as shown in Fig. 3.11. The pair-correlation functions are shown only for the interchain monomer pairs, i.e., all monomer pairs that do not belong to the same chain. It is important to point out that strictly speaking, we only need to exclude interactions of all monomers of a chain within the persistence length. However, if we consider the interaction of the same chain monomers smaller than the persistence length, there could be additional structures in the profile associated with the preferred distances along the polymer backbone which is caused by a combination of an average bond distance and the angles of the bending potential.

In the melt the g (r) has a liquid-like profile of regularly spaced broad peaks that correspond to neighboring shells of particles. Upon cooling below the crystallization temperature, these peaks become more pronounced and bi-modal, which is generally associated with more positional order. The g (r) does not give a clear picture due to the fact that it is averaged over all orientations. To examine the structure in detail, a parallel and perpendicular correlation functions denoted as g ⊥ (r) and g (r), respectively, are computed.

To compute these correlation functions, we assign to each monomer a local orienta- tion, which is the direction vector between the position of two bonded neighbors along the chain. For the end monomers we use the bond direction. For the perpendicular paircorrelation function g ⊥ (r) we compute a two dimensional distribution by considering a pair of monomers (i, j) when the projected distance along the direction ûi is less than σ and r is the distance perpendicular to ûi . In other words g ⊥ (r) is the correlation of particles confined between two parallel planes separated by 2σ, with the position and the orientation of the slab depending on one of the monomers.

A similar approach is used to compute the parallel pair-correlation function g (r) where r is the projected pair distance along the direction ûi . This calculation is limited to pairs of monomers that have a perpendicular distance less than σ. This is a scenario in which we consider a one-dimensional correlation, where the volume of a cylinder is considered along the unit direction of a monomer i. To compute both these pair-correlation functions the pairs (i, j) and ( j, i) need to be treated separately. Both these correlation functions are normalized to reach unity in the limit of large distances.

The peaks and crests observed in g (r) are more pronounced in g ⊥ (r) as shown in Fig. 3.11b. The perpendicular correlation distribution suggests that monomers locally have a hexagonal arrangement and hence, a similar type of ordering can also be expected in the aligned chain segments of the semi-crystalline state. The g (r) function presented in Fig. 3.11c reveals the correlation between the neighboring chains. The locations of the minima in the pair functions are correlated by integer multiples of the bond distance. It should be highlighted that the distributions shown have contributions from amorphous, crystalline and interface domains. In the melt, there are no obvious correlations to identify as the profile apart from a typical liquid-like behavior. However, at low temperatures the profiles are compatible with aligned, hexagonal arrangement of polymers with staggered positions.

A different way to characterize the local structure at the monomer scale is the local orientational order. A common approach to perform such an analysis is to divide the simulation box in cells. Each cell can contain several bonds and the local bond orientation within each cell is computed using the nematic order tensor Q αβ . The strength of the local orientational order is given by the largest eigenvalue S obtained from the

Q αβ Q αβ = 1 N i 3 2 b i α b i β - 1 2 δ αβ (3.4)
where b i is the unit bond vector, α, β are associated to the three directions x, y and z, and δ αβ is the Kronecker delta [124]. The nematic tensor is normalized by the number of bonds in the cell. We take a similar approach, but instead of cells we define a local spherical environment for every bond in the system. The position of the bond is considered to be at the mid-way point between the two corresponding bonded monomers and the radius of sphere is chosen to be 2σ. The main purpose of taking this approach is that the volume of the sphere remains unchanged, whereas in the approach where the simulation box is divided into cells, the shape and the size of the cells would depend on the simulation box.

The local orientational order can be used to define the degree of crystallinity χ of the system as the fraction of bonds in the system that is identified as locally crystalline, S ≥ 0.8 [131]. The results of crystallinity as function of the imposed temperature for both long and short polymers at different cooling rates is shown in Fig. 3.12. For the fast quenching rate the results are very similar for both the chain lengths. However, during the slow cooling nearly 60% of crystallinity is achieved for short chain systems, whereas in the case of long chains at the same rate of cooling only close to 50% is obtained. This reduced value in long chains can be associated with the slow dynamics.

We use the information of the local orientational order to perform the cluster analysis in a method similar to Morthomas et al [142,143]. To identify crystallites we first determine the crystalline candidates that are characterized as bonds with a local orientational order S ≥ 0.7. We choose a bit weaker value of S than that we use for χ. On these selected bonds, we apply a cluster algorithm in which the bonds belong to the same cluster if they are within a distance of 2σ and the relative angle between the bond directors is less than 3 • . Finally, only if a cluster has more than 10 bonds it is considered to be a crystallite.

The evolution of the number of crystallites and the average size of the crystallites as a function of the imposed temperature is shown in Fig. 3.13 for both the chain lengths. As the temperature reduces and reaches the crystallization temperature, both the above mentioned measurement increases rapidly. Thereafter, the growth-rate as function of the imposed temperature of the average crystallite size becomes more or less constant for both the chain lengths. The growth of crystallites is through the addition of small chain segments and monomers. During the slow cooling crystallites merge which explains the decrease in the number of crystallites below the reduced temperature of 0.75. For the fast quenching rates, the results are similar for both the chain lengths, we see a difference only for the slow cooling rate in the number of crystallites formed. In the long chain system there are less crystallites per unit volume compared to the short polymers, which was along the expected lines based on the χ measurement. Snapshots of the melt and the semi-crystalline configuration for a short chain system are shown in Fig. 3.14. The color coding used in this image is based on the local bond orientation order S. The melt configuration Fig. 3.14a is subjected to quenching at a constant rate 10 -6 τ -1 to obtain the semi-crystal Fig. 3.14b. The ordered regions formed in the semi-crystal are randomly oriented and have an average of nearly 200 bonds per crystallite.

Strain Induced Crystallization

Having analyzed the structure of both short and long chain systems during cooling, we apply the same type of measurements to investigate the structure of both the systems under strain. To study the strain induced behavior of the systems, we perform uni-axial deformation at different constant true-strain rates ranging between 10 -3 and 10 -6 τ -1 , close to the crystallization temperature. During the uni-axial deformation at a constant temperature the system is stretched along one direction, while pressure in the lateral dimensions is independently maintained at the same atmospheric pressure as before. Note that, similar to the slow cooling rate, the slow straining rate considered here is also unrealistically high for real experimental conditions. The configurations which are subjected to uni-axial de-formation are obtained from the slow cooling-rate ( Ṫ = 10 -6 τ -1 ) simulations at different temperatures near the crystallization temperature. The configurations obtained are directly subjected to straining, as the equilibration time for long polymers are beyond our computational resources and to make the comparative study of both chain length easier. The results presented in this section are based on single realizations, but we have verified that these results are representative using independent realizations. A series of snapshots at different strains is shown in Fig. 3.15 for a long chain system at a constant temperature of T = 0.83 under a slow true-strain rate of 10 -6 τ -1 . Conformations of an individual polymer under the same straining rate are presented in Fig. 3.16 at various strain values. This sequence of snapshots suggests that the initial stage of straining results in a weak deformation and rotation of polymer bonds. Upon further straining, the short straight segments stretch and continue to re-orient along the straining direction. Finally, the straight segments elongate along with the removal of folds and other defects. This observation is at odds with the findings of Yokota and Kawakatsu [47] for a single chain model of isothermal nucleation of polymers. The stress-strain plots for both chain lengths at different temperatures for two constant true strain-rates are shown in Fig. 3.17. The measured stress is averaged over short timeintervals and the stress-strain data is fitted with an unrestricted continuous piece-wise linear function with at most three segments to guide the eyes. For comparison, the result at the temperature T = 0.76 for the fast strain-rate is included in the same plot, which is already in the semi-crystalline state before straining. The remaining configurations obtained from slow cooling simulations are either in a melt state or weakly ordered as they are not subjected to an extended equilibration. The response of the stress as function of the imposed strain shows three stages. In the first regime, the short polymers have a very weak linear response that is more noticeable for the long chains. In the second regime there is a steep increase of building up the stress, followed by a steady increase in the third regime. Whether or not those regimes can easily be identified depends on the polymer length, the strain-rate and the temperature. Fig. 3.17a corresponds to the short polymers with the strain-rate of 10 -5 τ -1 . A weak linear response of the stress with respect to the imposed strain forms the first regime, and this corresponds to a melt system. Next, the disordered or weakly ordered polymers changes to more ordered structure in the second regime where the stress increases rapidly. In the third regime, the strain response of a semi-crystalline system can be seen. At high temperatures, both fast and slow strain-rates, 10 -5 τ -1 and 10 -6 τ -1 , respectively, show only the first regime, indicating that the system remains in the melt phase. The first regime is not present for the low temperature systems at the slow strain-rate as shown in Fig. 3.17b. These systems are already in the semi-crystalline phase as the polymers have more time to relax and adjust compared to the fast strain-rate. Hence, we observe a strain assisted transition. At T = 0.82 the system does not crystallize within the simulation time. At T = 0.83, the dependence on the strain rate can be seen for the short chain systems.

The results of the long chain systems are presented in Fig. 3.17c,d, for fast and slow strain-rates, respectively. The same observations regarding the temperature and strain-rate dependence can be made for long polymers as well, but the first regime is more pronounced and applies to systems well above the crystallization temperature. These results are quantitatively in agreement with other simulation [16,144] and experimental [145] observations of the stress-strain plots at different temperatures.

We extract the strain at which the stress response changes from the first regime to the second and refer it as "critical" strain. This strain value for various intermediate temperatures for both the chain lengths and the strain-rates is shown in Fig. 3.18 and a similar type of observation can also be seen experimentally [118]. For the systems at low temperature, the stress response is immediate since the systems are partially ordered and hence the transition strain would be a non-zero value. For the short chains at high temperature, the system remains in the melt and critical strain does not exist. For the long polymers the critical strain value is large, and for the high temperature system this strain might exist outside our simulated strain-range. In Fig. 3.18 we can observe that, the strain values for the fast strain-rates are always higher compared to the slow rates. The larger strain values found at faster strain rates appear to be counter intuitive. However, the actual time it takes to achieve crystallization is shorter. In other words, a faster crystallization is observed for faster strain-rates. It should be noted that, the initial configurations that are subjected to uni-axial deformation are not relaxed and hence, a long equilibration will affect the range of temperatures for which critical strain is obtained. To illustrate the effect of equilibration, we consider the short polymers, because the time to fully equilibrate the long chains is beyond our computational resource. The initial configurations, which were subjected to straining, are examined in an isothermal simulation at the same pressure as mentioned before. The result of crystallinity as function of simulation time is shown in Fig. 3 The change in structure can be characterized by pair-correlation functions excluding the intra-chain monomers, similar to the case of cooling. The results for the slow straining rate at the temperature T = 0.86 for the long chains are shown in Fig. 3.20. In these results, we use the strain direction as the reference, for parallel and perpendicular correlation functions. The results are similar in shape and structure to those observed during cooling, although in the cooling case we had considered the local chain direction as the reference. However, only at close proximity, the distributions g ⊥ (r) and g (r) have minor differences that stem from the nearest neighbors. Overall, the same observations can be made with respect to the local structure that points at aligned, hexagonally arranged polymers with staggered positions. The equivalence of the local structure between cooling and straining, does not extend to large length scales. During cooling the system will form various crystalline regions with random orientations, whereas under strain the ordered regions will align along the direction of strain. In Fig. 3.20 at the strain value of 1.2, the intermediate structure between melt-like and semi-crystalline behavior can be observed. This falls in the second regime, which is a sharp increase of the stress as function of imposed strain.

We measure the crystallinity during straining by using the same definition as above section 3.3. The results are presented in Fig. 3.21(a, c), for slow strain-rate for short and long polymers at different temperatures. The amount of crystallinity for both chain systems during straining exceeds that of cooling. Hence, the applied strain not only facilitates stretching but also enhances the amount of crystallinity.

A different local order is computed is the second Legendre polynomial P 2 [146] of the cosine of the angle between the unit-direction of a bond and the strain direction

P 2 = 3 2 b • ε - 1 2 , ( 3.5)
where b is the unit-direction of a bond and ε is the straining direction. The results are presented in Fig. 3.21(b, d), for slow strain-rate. We can observe that for the long chain systems the initial response of the stress, i.e., the first regime, mainly corresponds to a reorientation of the bonds before the actual crystallization begins. However, for the short chain systems under the slow strain-rate the first regime is not present as shown in the stressstrain plot of Fig. 3.17b and the same observation can be made for the P 2 measurement. Similar to the crystallization during slow cooling, the stems can be identified during straining. The results of slow strain-rate for short and long chain lengths are depicted in Fig. 3. 22(a,c). At low temperatures (T ≤ 0.82), some very short stems are already present in long chain systems. These short stems will grow upon straining, which induces crystallization. After reaching an average stem length between 10 to 15 units, about 50% of all the particles are contained in a stem as shown in Fig. 3. 22(b, d). The rate at which the stem length grows decreases, but the average length continues to increase steadily. At the largest strain, nearly 80% of all the particles are part of stems. The same behavior is also found at the higher temperatures, with the exception that at small strain where the transition from first to second regime occurs. For the highest temperature T = 0.96 that we have considered, a large number of short stems are formed with less than 20% of the particles in stems with an average length of less than 5 units. This is in line with our observation that the system remains in the melt-like state according to all the above measurements that we have considered.

The behavior of short polymers characterized by stem length and fraction particles in the stems are similar to long polymers at high temperatures (T > 0.83). The short polymer systems remain in a melt-like state with the average stem length of less than 5 units. However, the lower temperature systems only show two regimes, the second and the third. Here, the average stem length and the fraction of particles within stems reach higher values as compared to those of the long chain systems. At maximum extension the stem length reaches about 45 to 50 units and about 90% of all the particles are part of stems. The shorter polymers have fewer entanglements by which the dynamics and the structure formation is hindered. This leads to more noticeable structure, like the higher average stem length and the crystallinity. For the systems at higher temperature lack of entanglements allow for faster dynamics and the system remains in a melt state.

The strain induces the formation of different sizes of crystallites. A chain can connect to the same crystallite by looping, such a chain conformation formed is generally labeled as a hair-pin structure. In this study, we define such a hair-pin formed by two straight segments that belong to the same chain and are found in the same crystallite, but with opposite directions. The results for both chain lengths at the slow rate of deformation are shown in Fig. 3.23. Since there are fewer crystallites in the system, we observe initially more number of hair-pins. Upon straining the number of crystallites increases and the number of hairpins per polymer decreases. We can understand the formation of hair-pins in connection to other measurements such as stem length and crystallinity as function of the imposed strain.

In the regime where the average crystallinity and the stem length are very small, there are no crystallites formed so no hair-pins are found. At larger strains, when the crystallinity and the average stem length both increase and the number of stems tends to decrease the hair-pins dissipate slowly. This is an expected result, as this measurement is proportional to the chain length with the assumption of the fixed entanglement length.

When we compare the number of hair-pins per polymer formed during slow cooling with that of straining, we can see that the number of hair-pins increases with the decrease in temperature, due to merging of clusters as shown in Fig. 3.24. The merging of clusters can be inferred from the plots of the average number of crystallites and its average size as function of the imposed temperature. During the slow cooling rate, at low temperatures, both short and long chain systems show the decrease in the number of clusters, whereas the for long polymers at different temperatures.

average size of a cluster increases continuously.

We examine the structural changes at the polymeric scale by the means of end-to-end distance of monomers and the size/shape of the polymer. The size/shape of a polymer can be characterized by the radius of gyration tensor

R g,αβ = 1 N N i=1 r i α -r C M α r i β -r C M β , (3.6)
where i is the position of a monomer with respect to the center of mass r C M . Since we are deforming our simulation box along one direction, our system is elongated rectangular box. Therefore, we will measure the radius of gyration components that are parallel R g, = R g,zz and perpendicular R g,⊥ = R g,x x + R g, y y /2, to the straining directions. The results of both the end-to-end and the radius of gyration are presented in the Fig. 3.25 as a function of imposed strain. As expected, the chains stretch along the straining direction and reduce along the perpendicular direction. For short chains, the average end-to-end distance is about 40σ, which is around 80% of the maximum chain elongation that is possible. This chain extension is reflected in the values of the radius of gyration, which shows the increase of anisotropy of the overall polymer shape aligned with the straining direction. The relative anisotropy as a function of the imposed strain for both the chain lengths during the slowest strain-rate is shown in Fig. 3.26. For the short chain system at high temperatures T ≥ 0.83, this measurement does not change with strain, similar to the stress-strain curves. At these temperatures, the system remains in a melt-like state and the mobility is large enough to eliminate any induced deformation in the overall polymer shape. In the case of long polymers the shape and size change significantly due to the imposed strain at all the temperatures for slow straining. The time scale for which we have performed the simulations are too short for the long chains to disentangle even at the highest temperature considered here. Hence, the entangled polymers are forced to respond to the imposed deformation. Only at the highest temperature, T = 0.96, we see a different behavior of the stress response with respect to the imposed strain. However, it is inconclusive whether the system at this highest temperature also forms a semi-crystalline phase at even higher strains or that it remains in a non-crystalline state with some weak local orientational order.

The overall picture that emerges for the strain induced ordering for long polymers at temperatures above the crystallization point is that, the reorientation and ordering of polymers occurs in three stages. The first stage corresponds to the reorientation of the bonds upon stretching of the individual polymers. This takes relatively little effort as indicated by a weak linear stress-strain behavior. In the second stage, the fast increase of the stress, accompanied by the increase in crystallinity and straight segments in the chain is due to the entanglements present in the system that restricts the simple reorientation of bonds found in the first stage. The length of the straight segments or stems formed in this model is about 10-20 monomers. Hence, the continued straining of the polymers pulls the polymers locally taut. At the time the reduced mobility of monomers allows for the formation of crystalline regions, the hair-pins structures start to appear in both long and short chain systems. This relatively short domain ends according to our definition when the crystallinity is about 50%. Upon further straining, the stress continues to increase gradually, but with a smaller gradient compared to that of the second stage. The same is the case for various other measurements. The fact that the number of hair-pins decreases shows the re-organization that is taking place within polymers. This also suggests that the number of entanglements in the system is also reducing. For the systems below the temperature T ≤ 0.82, the same observations can be made, expect that the first regime is very short or does not exist. As mentioned earlier, the results presented in this section, are based on single realizations only. We have performed the same straining simulations in different directions, using the same configurations in the case of long chain systems, or used different realizations in the case of short chain systems. This results in a minor shift of behavior to either smaller or larger strain values. This a finite size effect, but the overall behavior remains the same. It is important to note that the results would be affected by longer equilibration time before applying uni-axial deformation. By surpassing this equilibration, we can focus on the response of a melt-like system for the short chains as well and confirm the correctness of the overall behavior that has been illustrated.

Conclusions

We have simulated a coarse-grained model inspired by PVA, where we focused on temperature induced crystallization via cooling from the melts and strain induced crystallization upon uni-axial deformation of the system, for temperatures at and above the crystallization temperature. For both the studies we have considered polymer systems with chain lengths of 100 and 1000 monomers.

For the case of temperature induced crystallization, different cooling rates were used and the model exhibits either a glass transition during fast quenching or a transition to a semi-crystalline state under slow quenching. Except for small differences in the density, the transition temperature and structural properties do not change for fast quenches. However, for the slowest quench rate considered in this study, the effect of entanglements in the long chains lead to slow dynamics which in-turn results in lower crystallinity and fewer crystallites per monomer. Here, the average stem length and crystallite size remain similar for both the polymers.

In the case of the long polymers, the self-interaction over the periodic boundaries cannot be excluded. However, if such interactions occur, they are separated by several hundreds of monomers along the polymer backbone. Hence, it is not expected to have a significant impact on the crystallization process.

During the uni-axial straining similar arguments can be made. With the given entanglement length, multiple entanglements with different polymers would screen a direct interaction of a possible self-entanglement. Upon inspecting various trajectories, we did not find any such problem.

The strain induced crystallization observed for both the chain lengths, occurs in three stages. In the first stage, polymers are still flexible and mobile thereby allowing for the local reorientation of polymer bonds. In the second stage, rapid crystallization was seen and at the last stage, the local structures, such as hair-pins, were removed from the system obtaining almost fully elongated chains. If we compare our model with the single and multichain models [47], we find that aspects of both these models, i.e., removal of loops and other defects within the individual polymers and the lateral growth of crystallites were observed.

It is important to note that the strain-rates we have applied would be very large under experimental conditions. This is due to the limitations available in the CPU-resources. A slower strain-rate, will obviously reduce the temperature range at which the strain induced crystallization can be observed. However, the qualitative picture is expected to remain.

In this particular model and for the strain rates employed, the short chains stretch to nearly 80% of the maximum extension at the largest strain. The trend, however, suggests that this will continue. For long polymers this is unclear, as the end-to-end distance only reaches upto 25% of the maximum extension in our simulations. Whether full elongation can be achieved for long chains cannot been concluded.

Chapter 4

Coarse-grained simulations of cis-polyisoprene

Introduction

Isoprene occurs naturally in various plant and animal materials, like natural rubber, terpenes, and camphor. In second half of the 19th century, isoprene was discovered along with its correct structural formula [147]. Natural rubber (NR) is mainly composed of cispolyisoprene, whereas isoprene rubber or synthetic rubber (IR) can contain primarily cis-PI or a blend of cis, trans and vinyl units depending on the synthesis. NR is extracted from the sap of the Hevea Brasiliensis tree and IR can be produced by polymerization of 2-methyl-1,3butadiene with a Ziegler-Natta catalyst [148]. However, there are several other methods to synthesize IR, like following the method by Rauline [149] used in the Michelin Tire Company. Both natural and synthetic rubbers are often mixed with other materials, like vulcanizing agents, oils, and nanofillers to enhance mechanical and physico-chemical properties [150].

The chemical composition of both IR and NR are very similar. The results from the literature show that in both vulcanized and unvulcanized forms NR has a higher tensile strength and tear resistance compared to synthetic analogues such as IR, polybutadiene (PB) and styrene-butadiene rubbers [151,152]. There are several studies investigating the SIC in rubber and various explanations have been proposed to describe the superior performance of NR. Among them is the presence of impurities in NR, such as proteins, phospholipids, carbohydrates and fatty acids [153]. These impurities in NR can act as permanent entanglement nodes which in turn enhance crystallization [154]. Candau et al. [22] and Toki et al. [28] examined the stress-strain relationship between NR and IR at various temperatures. The results from these studies suggest that vulcanized IR crystallizes during straining at lower temperatures and at larger stretching ratios compared to those of NR. Note that the SIC in unvulcanized IR occurs only at very low temperatures. Tanaka [153] examined the effect of non-rubber substances on thermal induced crystallization and found that immiscible stearic acid can be used as a nucleating agent in NR, while Kohjiya et al. [155] demonstrated that the addition of stearic acid had no influence on SIC in both NR and IR at room temperature. Many authors refer to micro-structural defects and a lower stereoregularity in IR to explain the lower crystallinity.

Polyisoprene will crystallize only if the polymer forms long sequences with stereoregularity [147]. These steroregular sequences must be linear structures composed of 1, 4-, 1, 2or 3, 4-isoperene units. If the chain is composed of 1, 4-, they must be either all cis or all trans. If the chain is composed of 1, 2-or 3, 4-isoperene units, they must be either syndiotactic or isotactic. Only cis-and trans-1, 4-polyisoperene can exhibit crystallinity, whereas 3, 4-isoperene is reported to be atactic because it is amorphous.

Computer simulations and modeling enable us to provide an explanation for experimental observations of the macromolecular structure and dynamics that originate at the microscopic level. They form an important tool to study polymers for physicists, chemists and engineers. To better understand the microscopic mechanisms of the experimentally observed macroscopic phenomena, computer simulations can be useful due to its ability to perform a qualitative analysis at the molecular scale. However, there are many challenges associated with simulation studies. For example, the accuracy of the force-fields for which the parameters are mainly obtained from ab initio quantum chemical calculations. Such simulations are performed at specific thermodynamic conditions, for instance at a constant temperature and pressure representing a melt-phase of small molecules like dimers, trimers or oligomers. Therefore technically the simulations must be performed using those specific conditions only. However, in practice the atomic force fields are not treated as temperature dependent as it is assumed that the potentials do not change too dramatically. Consequently, extensive testing and validation will be required if such force-fields are used for polymer simulations. Another problem is related to the time and length scale as briefly discussed in Chapter 1. Despite these challenges, computer simulations remain a very powerful tool due to its ability to study the microstructures in the existing polymer and also to predict new microstructures that can be useful in the development of new polymers. They can help by providing meaningful explanations for mechanisms involved in various processes, like the crystallization as mentioned in the introduction of Chapter 3.

Polyisoprene systems have been investigated by various research groups using different atomic, UA and CG models. Dodge and Mattice [156] used a modified CHARMM force field [157] to study trans-PI confined in crystalline channels and observed that conformational changes occur along the polymer backbone. Moe and Ediger [158] studied the local dynamics of cis-PI melts using the CFF91 potential. To understand the dynamics of the methyl group in polyisoprene at low temperature, Alvarez et al. [159] performed MD simulations using a Polymer Consortium Force Field. They found a good agreement between the static structure factor computed from simulations and that of neutron scattering experiments using spin polarization and also calculated the density of vibrational states of the methyl groups. Faller et al. [START_REF] Faller | [END_REF] developed a force field to study mono-and polydisperse oligomer melts of trans-1,4 PI by means of a combination of ab initio quantum chemistry and an automatic generation of empirical parameters. This atomic model was validated by comparing it with scattering experiments. Harmandaris et al. [160] refined the united-atom model developed by Smith and Paul [161] for cis-polybutadiene to study the density for different chain lengths as a function of temperature and self-diffusion properties in cis-PI. All these force fields have been used to study various local structures and bulk properties for nearly two decades.

Recently, Sharma et al. [150] performed MD simulations of both PI and PB systems by modifying the existing atomic and UA models introduced by Faller et al. [START_REF] Faller | [END_REF] and Harmandaris et al. [160]. In this study Sharma et al. used an equilibration criteria based on the convergence of the potential energy of the system and also computed various static and dynamic properties, such as the end-to-end distance, radius of gyration, bulk modulus, density and thermal expansion coefficient. This study also includes step wise cooling simulations to compute the glass transition temperature.

Polymers have long relaxation times. Therefore, investigating phenomena at mesoscopic length scales using atomic or united-atom models requires more computational resources than that are currently available. This is due to the number of degrees of freedom involved in the atomic and united-atom simulations. In the case of atomic simulations, each atom has six degrees of freedom which correspond to three spacial coordinates and three velocity components. Coarse-grained models have been used to perform simulations at the mesoscopic scale by retaining the most relevant degrees of freedom to investigate a particular phenomena. The idea of coarse-graining has been in use for quite some time. Here a group of atoms or united-atoms are mapped onto a single CG bead. Since coarse-graining by definition reduces the number of particles in a system, the number of degrees of freedom is also reduced and hence, bigger length and time scales can be accessed.

Li et al. [START_REF] Li | [END_REF]162] performed atomic simulations using COMPASS [163] and obtained the CG potentials using the IBM 2.3.2. The geometrical mapping from the full atomic to the CG scale is done by selecting the center of the single bond that connects two monomers as the mapping point of the CG bead. This mapping scheme leaves the ends of the chain with a very anisotropic environment. The CG bond potentials were modeled using a harmonic potential, while the angle and the non-bonded potentials were tabulated. In the same study, cross-links were modeled by introducing two additional beads, the cross-link site and the cross-linker bead, where the former can have either 2 or 3 bonds and the later can have a maximum of 2 bonds, respectively. The cross-linking process was done by a CGMD simulation based on a distance constraint between the cross-link site and the cross-linker bead. Further, they developed a hierarchical multi-scale model where the information from the atomic simulations is passed to a continuum constitutive model to predict the mechanical properties. The main idea in this hierarchical scheme is to decompose the material into 2 structures, a cross-link network and a free-chain network representing hyperelastic and viscous models. The cross-link network is characterized by the number of cross-links and the length between the cross-links. A well established Arruda-Boyce model [164] was used to describe the elastic behavior of the hyperelastic network. The viscous properties of polymers are considered to be a dynamic property of free-chain network. They performed straining simulations using the finite element analysis.

Ohkuma and Kremer [165] performed united-atom simulations to obtain the CG potentials using the IBM. In this study they compared angle and bond distributions for different CG mapping schemes. They selected a mapping scheme with two beads per monomer as they found both the distributions nearly independent to each other. Using this mapping scheme they obtained two CG-potentials, with and without pressure corrections and computed various properties like density and entanglement length.

More recently Shahidi et al. [101] performed MD simulations with a UA model at different constant temperatures and developed CG potentials using the inverse Monte Carlo method 2.3.3. In addition, they also developed another CG potential which included an order parameter that is based on the nearest-neighbor distribution to capture the local dynamics. It was shown by Henderson [166] for the mono-component case and generalized later by Rudzinski and Noid [167] for the multi-component and intermolecular interactions that the relationship between pair potential and RDF results in a unique solution. However, this relationship may not always be unique in the IBM. Soper [88] has shown that different potentials may in some cases result in similar RDFs. Because of this, both methods the inverse Monte Carlo and the iterative-Boltzmann, can produce different results while having very similar RDFs. The choice of selecting the mapping center of the CG bead is not unique, there can be several options as shown by Ohkuma and Kremer [165]. In Fig. 4.1 different CG schemes for polyisoperene system are depicted.

The most intuitive way to define the center of a CG bead is the center of mass of all the atoms in a monomer as shown in Fig. 4.1a. This mapping scheme is used by Shahidi et al. [101] and the probability distribution of the bond lengths results in a double peak, as the carbon-carbon single bonds can easily flip from one torsional state to the other. Ohkuma and Kremer [165], and Li et al. [81] demonstrated that the Type-1A mapping option results in a correlation between the bond length and the angle, which are therefore not independent.

The Type-2A mapping method is shown in Fig. 4.1b. All the atoms between carboncarbon double bonds are considered and the center of the CG bead is located at the COM of these atoms. This mapping scheme has not been used in literature as the bond length distribution does not result in a clear single peak due the effect of the torsional rotation around the carbon-carbon single bond. We have considered this mapping scheme as it is a rough approximation of the Type-2B method which is widely used in literature.

The mapping scheme Type-2B is shown in Fig. 4.1c. This mapping method is adopted by Faller et al. [START_REF] Faller | [END_REF] and Li et al. [START_REF] Li | [END_REF]162] to avoid the double peak distribution of the bond lengths. The rationale is that, when the bond length distribution has a clear single peak, it can easily be represented by a harmonic potential. Also the position of the CG beads do not depend on the dihedral angle formed between the carbons of CH CH 2 CH 2 C . The angle distribution shows a double peaked structure. It is mentioned in the study by Ohkuma and Kremer [165] that the Type-2B mapping scheme can also have correlations between the bonded variables. This means strictly speaking that the bond-angle potential is required to model such interactions. However, the bond and angle potentials are generated independently in all the studies available in the literature using this mapping scheme. Therefore, in the CG model the bond and angle correlations do not exist. Hence, the measured structural distributions can be different in both CG and atomic simulations.

The Type-2C mapping option is depicted in Fig. 4.1d, where the location of the CG grain is at the COM of atoms CH 2 CH 2 . The center of the CG bead will be very close to that of the Type-2B option.

The Type-3A and Type-3B mapping options have two beads per monomer as represented in Fig. 4.1e and f, respectively. In Type-3A one bead constitutes CH 2 C CH 3 atoms with the mapping center located at the COM of these atoms and the second bead is located at the COM of the CH CH 2 atoms. The study [165] shows that the bonded variables, angle and bond length are correlated.

Ohkuma and Kremer [165] selected the Type-3B mapping model. The center-of-mass of two carbons in CH 2 CH 2 forms one bead and the COM of three carbons in CH 3 C CH the other bead. This mapping scheme was selected due to the fact that the angles and bond lengths were found to be almost independent of each other.

A CG model constructed by using the IBM assumes that the bonded variable sets like bond lengths, angles and dihedral angles are independent with respect to each other. Hence, it would be relevant to verify these variables are in fact independent. The study in Ref. [165] only checks the dependence of two bonded variables, bond lengths and angles. They observed that only the Type-3B scheme satisfies this criteria.

As per our knowledge the mapping schemes Type-2A, Type-2C and Type-3A have never been used. We can conclude from these studies that a good CG mapping depends on many things, like the method used to obtain the optimized CG potentials and the goal of the study. Also, it is worth noting that the iterative Boltzmann inversion and the inverse Monte-Carlo methods considers only the structural measurements like the radial distribution function and the dynamics of the system is not accounted for during the development of CG potentials. However, the validation procedure for any CG model remains the same, irrespective of the above mentioned factors, which is how well can the CG model reproduce the static and dynamic properties with respect to the reference MD measurements.

The goal of this study is to obtain a CG model to understand the strain-induced/assisted crystallization of cis-polyisoperene by using statistical trajectory matching method. To achieve this target we perform atomic simulations at a constant temperature. The MD configurations are stored regularly at a desired timestep. Then, different CG mapping schemes are employed to obtain respective reference CG trajectories. These reference trajectories are subjected to STM fitting process and optimized CG potentials are obtained for different CG mapping options. Later we perform DPD simulations using the CG potentials. In the studies mentioned in the literature survey above, the MD simulations using atomic or united-atom models are performed either using Gromacs [168] or a commercial package like Material Studio [163,169]. In the simulations performed using the Gromacs code, all the bonded atom types are considered as rigid bonds. In this work, we use a generalized force field called AMBER where rigid bonds are not used. These potential parameters are obtained mainly for small organic molecules using quantum chemistry calculations.

Atomistic Model

The total potential energy is the sum of bonded and non-bonded interaction potentials. Bonded interactions include harmonic potentials for bond and angles, the proper dihedral interaction is implemented using the Fourier style and the improper dihedral interaction is not included in this work. The non-bonded interactions are modeled using the LJ potential with a cutoff radius of 9 Å and Coulomb interactions. The charges included in our model are not explicitly derived from the quantum calculations as described in the AMBER manual. We make a very crude approximation, that an united-atom entity is charge neutral, i.e., the sum of charges associated to a carbon atom and its attached hydrogen atoms is zero. Using this idea, we assign partial charges to every atom where all hydrogen atoms have +0.06 electron charge. Based on the number of hydrogen atoms attached to a carbon, a partial charge is assigned. The intramolecular non-bonded interactions for all the atoms separated by three or more than three bonds are considered. All 1,4 interactions are scaled down by the weight factor of 0.5 and 0.833, for LJ and Coulomb interactions, respectively. The force field parameters used in this work are provided in Appendix A.1. The timestep used in the atomistic simulations is 1 fs. The pressure and temperature of the system are controlled using a Nose-Hoover thermostat and barostat, respectively. We have applied periodic boundary conditions in three directions for all the systems.

Initial configurations

We generate atomic systems for 5 different chain lengths, 10, 20, 30, 40 and 80 monomers per chain using Medea software [170]. The initial configurations are generated with low density and have nearly 10,500 atoms each. After performing energy minimization on these initial configurations using a conjugate gradient method, it was observed that atomic systems had a variety of issues like • Presence of trans-monomers and • Different types of knots These often occur in real polymers, as the experimental systems are not always composed of 100% cis-monomers and they may also have knots. In fact knots appear naturally in long polymers [171] and play a vital role in crystallization. The energy required for a conformational change from trans to cis state was calculated from the isoprene spectra in 1960's. Subsequently, it is considered that the energy difference was not large enough to prevent cis to trans conversion [147].

Two examples of different types of knots found in the atomic systems are shown in Fig. 4.2. Some type of knots like in Fig. 4.2a can be easily un-knotted by stretching the polymer, whereas other type like in Fig. 4.2b can become tightly knotted upon pulling the ends of the chain. Although knots and trans-monomers do occur, their presence will lead to an unwanted second peak in the bond length distribution which can further affect the coarse-grained potentials obtained from STM process. Hence, it would be ideal to have a homogeneous (only cis-monomers) system without any knots. To achieve this homogeneous system, we take the following approach.

1. Energy minimization is performed by reducing the timestep and increasing the dihedral constant between the carbon atoms CH 2 CH C CH 2 in the polymer backbone 2. The trajectory obtained in step 1 is checked for the presence of any trans-monomers. If they exists step 1 is repeated with a reduced timestep.

3. A NVT simulation is performed using the same settings specified in the step 1 4. The NVT simulation is continued using the last configuration obtained in step 3, but the timestep is increased slowly to the desired value 5. After verifying that the trajectory obtained in step 4 has only cis-monomers, the pressure of the system is increased slowly until the density of the system reaches the desired average value in a NPT simulation 6. Then NPT simulation is performed at the desired timestep and pressure (one atmosphere) for further equilibration 7. Again the trajectory is verified to have only cis-monomers. If this is the case then the dihedral constant is slowly reduced to the desired value.

It was also observed that some chains had knots which were assisting the formation of trans-monomers. To avoid this problem, we identified the chains with knots and replaced the whole chain with a different one from another configuration. Thereafter the above mentioned process is followed to obtain a homogeneous system without knots. However, we have used a combination of bond length distribution, g(r) (excluding 1-2 and 1-3 particles) and visual approach to identify the knots in a coarse-grained configuration. To verify that our configurations are defect-free, we continued to perform the NPT simulation and employed the Type-2B mapping scheme on the MD trajectory. The coarse-grained atomic trajectory was used to measure the bond length, angle and g(r) distributions. In one case we found a longer tail in the bond length distribution as compared to the other atomic CG trajectories and also the g(r) had a small peak at a short distance. This meant that the configuration had some sort of defect which was not identified by the process we followed and hence, that system was rejected.

Upon verifying that the configurations did not have the above mentioned issues, they are subjected to longer equilibration. This simulation run was performed using 1 fs timestep, at a constant temperature and pressure, 350K and at 1 atmosphere, respectively. The chain lengths of 10, 20, 30, and 80 monomers are equilibrated for 80, 100, 150, and 300 ns, respectively. We have used only one realization for these systems. However, four realizations of chain length 40 monomers were subjected to each 150 ns of equilibration run.

Validation of the Atomistic Model

The above relaxed systems are used to measure various properties like the density, the radius of gyration and the end-to-end distance. Here we present a comparison between our results and literature findings.

The mass density measurement for the atomic systems of 5 different chain lengths at 350K are presented in Table 4.1. We can observe that the density of the systems increase as a function of chain length. This is an expected result, as the effect of end particles is higher for short chains. The difference in density for different chain length comes from the fact that the effective volume for the end particles of a chain forms a single-bond and therefore has a slightly higher effective volume compared to the internal particles that have two bonds. Similar observation regarding the polymer length scale can be done from the results of atomic cis-polybutadiene system [150] regarding the glass transition temperature measurement and also the fast field cycling NMR relaxometry experiments by Ref. [172].

To compare our model with other simulation models found in literature, the density of chain length 30 monomers is presented in Fig. 4.3. In the same plot experimental results are also included. The mass densities measured in our model underestimate the results found in the literature. The maximum deviation is less than 3 % when compared to experiment [173]. Similarly the density of our atomic system at temperature 298K is also lower compared to the atomic model used in Ref. [150], where they consider a system with the chain length of 32 monomers which is slightly longer than our atomic chain. We compare the two atomic models under the assumption that the chain lengths of 30 and 32 monomers are similar. At high temperature of 413K, our results are also lower than the UA model used by Shahidi et al. [101]. At temperature 350K the deviation of our measurement is less than 1 % compared to the work presented in Ref. [101].

Results are based on a single realization, except for the system of chain length 40 monomers, are used to compare our model with experimental data and other simulation results. Our model is within 3 % deviation when compared to the measurements found in the literature.

We have calculated structural properties of the atomic systems, like radius of gyration and end-to-end distance. The results for five chain lengths are presented in Table . 4.2. We can notice that the ratio 〈R 2 ee 〉 / 〈R 2 g 〉 ≈ 6, which for linear polymers is the characteristic value for an ideal chain. The end-to-end distance is computed between the end carbon atoms with three hydrogen atoms (cis-type) of the polymer backbone on both the sides of the chain. The end carbon atom with three hydrogen atoms on one side of the polymer chain can be clearly identified. But on the other side of the polymer chain, carbon in the methyl group can be found at two locations. One along the backbone of a monomer (cis-type) and the other is the side group (trans-type) attached to the single carbon atom which is part of the double bond.

Based on these results we can say that our atomic model is in good agreement with respect to the densities for the available literature data concerning simulation studies [101]. From the comparison between our measurements and the experimental data, we see that they are within 1 % deviation at high temperature. This result indicates that our model is in [173]. The measurements of a UA model used in Ref. [101] are shown as green squares and all-atom model from [150] is shown as blue square. The results from our model (AMBER potential) are the red triangles. The dotted lines are guide to the eye. good agreement with respect to the experimental result, given the fact that our chains have low molecular weight compared to the experimental systems and contains only one type of isomer.

Coarse-Grain models

As discussed earlier, the aim of this work is to obtain a CG model to study strain-induced/assisted crystallization in cis-polyisoperene using the statistical trajectory matching method 2.3.5. The production runs are performed at constant temperature and pressure. The purpose of these simulation runs is to obtain a MD trajectory where a configuration is stored at every 50 fs, which is the commonly used timestep (δt) in a DPD simulation that uses CG potential generated by STM [111]. The number of configurations required to obtain an optimized potential will depend on many factors like the chemical complexity of a polymer system, the cutoff distance used to obtain conservative and dissipative forces, and the selection of number of parameters to be optimized. For a simple system considered in [80], it is mentioned that the Bayesian optimization method converges the potential parameters pretty quickly. In our system, we observe that 10, 000 configurations are required to obtain a set of optimized parameters.

Our goal is to obtain a CG model for cis-polyisoperene that can be used to study the strain-induced/assisted crystallization behavior. To this end, we examined four mapping schemes, Type-1A, Type-2A, Type-2B and Type-2C as described in Sec. 4.1 and depicted in Fig. 4.1. The MD trajectory obtained from the production run of different chain length systems are applied to the different mapping options mentioned above, which we refer to as reference trajectories. To apply the STM method, appropriate interaction styles have to be selected.

Selection of interaction styles

As outlined in Chapter 2, there are several choices available to model the bonded interactions. The approach we take to select the appropriate style of bonded interactions is based on the shape of the measured structural distributions. The bond length and angle distributions are computed between two and three successive beads from analyzing the reference trajectories. The results of such calculations for different chain lengths using three mapping schemes are presented in Fig. 4.4. Note, the distributions for the Type-2C mapping are not presented here as they are very similar to the Type-2B scheme.

The bond length distribution has a bi-modal shape in the Type-1A mapping since two types of conformations are possible along the monomer backbone due to the changes caused by the torsional rotations. These conformations correspond to trans (φ ≈ ±180 • ) and gauche (φ ≈ ±60 • ) state formed by four carbon atoms of the monomer backbone. For the same mapping scheme the angle distribution has one peak with a skewed profile. To capture such behavior of bonded measurements we use quartic potentials to model the bond and angle interactions which are implemented in LAMMPS packages under the names of class2 and quartic interactions, respectively.

E bond (r) = K b2 (r -r 0 ) 2 + K b3 (r -r 0 ) 3 + K b4 (r -r 0 ) 4 , (4.1a) E ang l e (θ ) = K a2 (θ -θ 0 ) 2 + K a3 (θ -θ 0 ) 3 + K a4 (θ -θ 0 ) 4 , (4.1b)
where the r 0 , θ 0 are the equilibrium distance and angle, respectively. The K's are prefactors expressed in the units of energy/distance 2 for class2 bond style and energy for quartic angle style, respectively. All these eight parameters are optimized using STM to obtain a CG model for the Type-1A case.

In the case of Type-2A mapping, the bond length distribution exhibits a skewed shape profile. Similar to the Type-1A mapping, two types of conformations are possible even in the Type-2A approach due to the torsional rotation around the carbon-carbon single bond. The absence of two clear peaks in this case, as opposed to Type-1A, is the consequence of the location of CG bead which is the center-of-mass of all the atoms between the consecutive peak. The bi-modal angle distribution has the first peak at a similar location as in the case of Type-2A method, whereas the second peak is located at a larger angle.

It is important to note that in both Type-2A and 2B models the end atoms of the chains are neglected, which effectively means that we have one monomer less in each chain. This results in lower density as compared to the atomic system.

The conservative non-bonded interactions are modeled by a pairwise potential that depends only on the distance between non-bonded CG beads. The grains are implicitly assumed to be spherical. To avoid the non-bonded beads to come very close to each other different pairwise interaction styles can be used like, LJ, piece-wise polynomials and tabulated. More details on modeling the pair-interaction using piece-wise polynomials can be found in Ref. [110]. It is desirable to have more parameters to allow greater flexibility in the shape of the interaction potential. Tabulated potential is one such option where the forces are based on a linear interpolation between tabulated points [111]. One disadvantage of this type of interaction is that at small distance, there is a considerable uncertainty in the force because these distances occur very rarely in the reference trajectory [111]. Despite this drawback we use this tabulated method to model conservative non-bonded forces for all the above mentioned CG models with a cutoff of 20 Å. 4.3 shows the number of parameters used to obtain a CG potential using STM for various models. The pairwise conservative forces in all the models uses a tabulated style with 20 points where the force at first and last points are set to zero. The tabulated values are spaced with the mesh size of 1 Å. The models Type-2A and 2B use a tabulated angle interaction style with separation of about 5 • . We have also tested different number of points in the tabulated angle potential where the separation between each point is about 2.5 • . STM fitted angular forces obtained for these two cases were virtually identical.

For the non-conservative dissipative forces, we use the standard pairwise DPD expression as shown in 2.28 and 2.30 with a cut-off of 10 Å. The frictional parameter, γ is obtained analytically and remains identical for all CG bead pairs including the bonded neighbors.

The cutoff distances used for conservative and non-conservative forces are based on earlier studies using simpler models [80,110,111]. This means we assume that the settings used for relatively simple models are still valid and holds good for our model as well.

We have also used the Type-2C mapping scheme to obtain a reference trajectory and measure the bond length and angle distributions. These profiles have nearly the same shape as that of the Type-2B model. One difference that can be observed in the bond length distribution is that the profile is not as symmetric as in the Type-2B distribution. This difference is because of the location of the center of the CG bead in Type-2C mapping which includes the atoms in the carbon-carbon single bond along with its attached hydrogen atoms. These hydrogen atoms can vibrate very easily causing the discrepancy for the above mentioned behavior. However, the angle distributions for both the models are visually indistinguishable.

Coarse-Grained Potentials from STM

The reference trajectory obtained after applying the mapping schemes are used in the STM fitting procedure. We have one reference trajectory for every chain length that corresponds to one mapping. STM provides coarse-grained potentials which include the parameters for the bond, angle and pairwise interactions. Along with these parameters, a friction coefficient is also obtained. In Fig. 4.5 the angle and bond potentials for various chain lengths using different mapping schemes are shown. The potentials for the chain length of 40 monomers are obtained separately for four independent realizations and these potentials are combined to one unified potential.

The shape of the bond and angle potentials obtained from STM are very similar for different chain lengths for each CG model. However, in the quartic bond potential for the Type-2A model, we can observe a systematic chain length behavior where the equilibrium distance is slightly increased with respect to the chain lengths. In this model, we can observe that the tabulated angle potential of the system with the chain length 20 monomers, is deviating from rest of the chain lengths at small angles. This is due to the fact that such small angles occur very rarely in the reference trajectory, causing a considerable uncertainty in the force at these small angles.

In the case of Type-2B model, for the chain lengths of 30 and 40 monomers the angular energy is slightly shifted upwards with respect to the other chain lengths in the range between 70 • to 120 • . In the case of longer chain systems it is observed that the duration of simulation plays a role in having slightly different angular energies. For example, the potential obtained after a long simulation of 300 ns has slightly more angular energy when compared to that of the potential obtained from a shorter simulation of 150 ns. This means our system with the chain length of 80 monomers requires considerably long time to reach equilibrium.

The pair potential obtained from the STM method has a very weak attraction in all the models irrespective of the chain lengths as shown in Fig. 4.6. However, it is important to realize that for a given model the conservative forces remains the same for all the chain lengths. At short distances, the pairwise conservative forces have some uncertainty as mentioned in the section 4.3.1. Due to this, some deviations are observed between the chain lengths for the same model.

DPD simulations

The DPD simulations are performed using the CG potentials obtained from STM method with a timestep of 50fs. All the CG models are subjected to 100ns of equilibration run followed by 100ns of production run in the NPT ensemble. The temperature is kept constant at 350K and the pressure at 1 atmosphere using a Berendsen thermostat and barostat, respectively. The CG simulations are performed for configurations with about 800 monomers for each model. We perform the simulations on small systems as a test case to observe the behavior of the structural measurements of our CG potential. We measure several structural and dynamic quantities to check the validity of our CG models. These measurements are compared to their MD targets obtained beforehand from analyzing the coarse-grained reference trajectories.

The comparison of the measured mass densities obtained from DPD simulations and its corresponding reference MD values for different models for various chain lengths are shown in Fig. 4.7. The densities obtained from CG simulations for different chain lengths are in good agreement with respect to the reference values for Type-1A model. The mass densities for all the chain lengths during the DPD simulations are within 1 % deviation from the reference MD values for Type-1A case. Note, the experimental density at 350K is computed by performing linear interpolation of the values obtained in literature [173] at temperatures 298K and 413K. The data for the Type-2A and 2B models are also presented in the same plot. We recall that in these mapping schemes there is one monomer less per chain. In the Type-2B CG model the density of the system is overestimated by nearly 3 % for the chain lengths greater than 40 monomers in the DPD simulations when compared with the reference measurements. However, the densities for different chain lengths are underestimated in the Type-2A CG model with respect to the reference measurements, where the discrepancy between the DPD and the MD values reduce as function of the chain length. The mass density for the short chain system has a greater deviation between the DPD and the reference MD values in the Type-2A model. For the chain lengths greater than 20 monomers, the density measured during the DPD simulations is within 4 % deviation with respect to the reference values.

The CG interactions obtained from the STM are different for the Type-2A and 2B variants, and this could lead to different persistence length in the polymers. The combination of these factors like persistence length, interactions (bond, angle and pair-wise) will affect the measurements of the CG system.

Structural properties

The structural measurements like bond distance and angle distributions are presented along with the reference CG MD curves in Fig. 4.8. The distributions shown are measured during the production run of the DPD simulations for all five chain lengths and the three different CG models. The reference MD distribution is the average of four independent realizations of the system with chain length of 40 monomers. This reference MD curve is selected as a representative for all the other chain length distributions, as they are very similar to each other.

The peaks of the bond length distributions are located at nearly the same distance for all the CG models irrespective of the chain lengths. The bond length profiles obtained from In Type-2A CG model, the bond length and angle distributions are well reproduced, but not all the peaks in the radial distribution functions are reproduced as shown in Fig. 4.9. This suggests that the Type-2A model cannot serve the purpose of our goal to study the crystallization in cis-PI inspired CG system.

Based on the structural measurements like bond length, angle and g(r) distributions for different CG models, we can say that the CG potentials obtained from STM is capable of reproducing the structures better in the Type-2B model than in the Type-1A and 2A. Since the goal is to obtain a CG model based on cis-PI to study the strain-induced/assisted crystallization, we have considered Type-2B model for which we have five CG potentials which correspond to five different chain lengths.

When we consider different models obtained for different chain length for the Type-2B case, we can see that CG simulations using 30 and 40 monomer systems are better than the other chain systems with respect to density and other structural measurements. The density for long chain systems (≥ 30 monomers chain lengths) is about 2 % deviating with respect to the reference measurement. The local structure measured from the RDFs is in accordance with the reference distribution. Based on these findings we can conclude that in the Type-2B model the chain lengths of 30 and 40 monomers are similar to each other. 10 and 20 chain length models suffer from end monomer effect, whereas the long chain system with 80 monomers suffer from slow dynamics. Therefore, CG models of 30 and 40 monomers are equally good in reproducing the reference measurements when compared to 10, 20, and 80 monomers for Type-2B case. We select the CG potential of 40 monomers in our model to build bigger and longer chain system. Further we measure the radius of gyration and end-to-end distributions for the small systems of different chain lengths using the three CG models in Fig. 4.10. The profiles for different CG models are similar for a given chain length. Flory's theorem states that long polymers in the melt are nearly ideal chains. In such case, it can be shown that the probability distribution function of the end-to-end distance is a Gaussian distribution. This does not depend on other structural measurements. We make such similar observations in our different CG models as well.

In order to compare the measured R g between the CG models and the reference MD measurements, we compute the average R g as shown in Fig. 4.11. It can be observed that the average measurements are systematically overestimated in the CG models in comparison with their reference values. In the MD case for longer chain systems the R g distributions suffer from poor statistics. This results in the poor estimation of the average values. Note that in the case of Type-2B model for the chain length of 40 monomers nearly 20 % deviation is found between DPD and MD measurements, where the MD value is obtained from four independent realizations.

We recall for the chain length of 40 monomers four independent configurations were each simulated and subjected to STM process. These four independent obtained CG potentials were combined to get one coarse-grained potential.

To overcome any finite size effects of the small system of 800 monomers (of the chain length 40 monomers), we generate a big system with 2500 chains by replicating the small system five times in x, y and z directions. The big system is simulated for 50ns and various structural measurements are computed. The bond length, angle and the radius of gyration distributions of small and big systems are visually identical to each other.

The probability distributions of angle and bond length obtained from the combined potential are used to generate long chain systems of 80, 160, 400, and 800 beads consisting of 100, 000 beads per system. The generated long chain configurations are subjected to an equilibration process where the repulsive forces are slowly introduced. Next, these systems are simulated at a constant temperature of 500K for 60ns, 200ns, 3000ns and 6000ns, respectively. Later, these equilibrated configurations are cooled down from 500K to 350K. After this, the systems are subjected to 50ns of isothermal and isobaric simulations. The densities of the big systems for different chain lengths are presented in Fig. 4.12. Note, the experimental density at 350K is computed by performing a linear interpolation of the values obtained from Ref. [173] at temperatures 298K and 413K.

The measured mass densities using DPD simulations for small and big systems of chain length 40 monomers are in agreement within 1 % discrepancy. The density of the big system increases with the chain length this is because of the effect of the end monomers that reduces in the long chain systems. However, after the chain length of 160 monomers the densities of the systems do not alter much. The mass density of the big system of chain length 40 monomers is nearly 3 % higher compared to the reference MD.

The probability distribution of the radius of gyration measured in the DPD simulation using the Type-2B combined potential are scaled with respect to the square root of the number of monomers as shown in Fig. 4.13. It can be seen that for the chain lengths below 160 monomers the R g profile shape is similar with a single peak. In the case of long chain systems of 400 and 800 monomers the R g distributions do not have a smooth profile as the simulation time required is much higher than what is considered in our study. The slight shift of the peak for the chain length of 40 monomers, indicates that these chains are still too short to show the full polymer scaling behavior for this measurement.

Dynamics

We have examined the ability of different CG models to provide an insight of the dynamics in the DPD simulations. The advantage of CG simulations is to access longer time and length scales compared to the atomic simulations. Here we present the mean square displacement (MSD) calculations of center-of-mass of the chains and the coarse-grained beads for different chain lengths of the three CG models in Fig. 4.14. The raw data presented here is for both DPD and MD measurements.

Note that the DPD simulations are performed using the CG potentials obtained for the respective chain lengths and CG mapping schemes. For instance, m010, m020 and so on corresponds to the chain lengths of 10, 20 and so on monomers. The CG potentials are obtained for each chain length after independently performing the atomic simulations and CG mappings.

We can see that the initial regime of the MSD curves are alike for both MD and DPD simulations. Beyond 1 ps of simulation time the DPD simulations are significantly different from the MD measurement. We can notice that for the MD measurements of the MSD curve, the COM of the polymers initially shows a ballistic regime followed by a regime where the polymers appear to be trapped. In the case of short chains it appears that the polymers escape the confined regime near 1 ns.

The COM of polymer in the DPD simulation moves more than in the MD simulation during the given time period of 1 ns, which indicates that the DPD simulations are much faster compared to the MD simulations within this time window (between ps and ns). A similar type of observation can also be made in the MSD for monomers.

Cooling Simulations

We have considered small systems (800 monomers) to perform cooling simulations. These small systems are plagued with various finite size effects and the measurements obtained will have big fluctuations. Despite this we perform cooling simulations using the CG potentials obtained for different chain lengths to get a quick idea about the behavior related to phase transition. For completeness, we present the results of CG simulations obtained from three CG mapping schemes (Type-1A, 2A and 2B).

We perform cooling simulations at different constant rates between 0.2 and 20 K/ns for the small systems from the temperature 350K to 150K for different chain lengths using various CG models at atmospheric pressure. The volume per monomer and average R g as function of the imposed temperature are shown in Fig. 4.15. It can be noticed that in the Type-1A model the volume per monomer and the 〈R g 〉 remain nearly the same as function of temperature for all the chain lengths.

In the Type-2A model the volume per monomer reduces with respect to the temperature with a constant slope, but the 〈R g 〉 does not change with reducing temperature. This suggests that the system does not undergo any phase transition in this temperature window. One possibility is that a phase transition might occur at even lower temperature than the temperature range presented here. Another reason for not showing a phase transition could be associated with the cooling rate. To test the hypothesis for some type of phase transition occurring at lower temperatures, i.e., below 150K we continue the cooling simulation from 150K to 1K at the same rate as shown in the plots which is 2K/ns and also to test the effect of cooling rate we perform slower cooling at the rate 0.2K/ns. We observe that the slope of volume per monomer as function of the imposed temperature remains constant until nearly 70K and below this temperature the slope becomes smaller. A similar observation can also be made for the slower rate of cooling where the slope change occurs at a slightly higher temperature around 90K. The 〈R g 〉 measurement fluctuates around an average value similar to Fig. 4.15d. At low temperatures (< 150K) the fluctuations for the average radius of gyration reduce. This means even if we consider the temperature at which the slope change occurs which is our estimate for the glass transition temperature, the Type-2A model underestimates the T g compared to the experimentally observed glass transition temperature of cis-PI which is about 200K [174].

In the Type-2B model, we see a phase transition for the chain lengths greater than 20 monomers which is indicated by the drop in temperature volume plot shown in Fig. 4.15e. The transition temperature for different chain lengths varies between 275K to 250K. The highest transition temperature is found in the 30 monomers chain length model followed by 40, 20, and 80 models. A plausible explanation for different transition temperatures can be found by looking at the angular potential in Fig. 4.5f. In the angular potential we can see that 30 model system has higher angular energy making it easier to move the angles from one local minima to another. This argument can be extended to 40 monomer chain length model as well. In the case of 20 and 80 models the angular energies are close at the first minima and the energy difference slightly increases at the second maxima, where the ,d,e) to the MSD of CG beads, respectively. Note that the CG simulations are performed using the CG potentials obtained for the respective chain lengths and the CG mapping schemes. The combined CG potential is used only for the chain length of 40 monomers system. 80 chain length model has higher angular energy. This means that in the 20 chain length model the angles can move easily to higher values as compared to the 80 chain system. In the Type-2B model, for 40 monomer chain length system we have considered four independent realizations, whereas for the rest of the chain lengths only one realization is considered. The transition temperature in the temperature-volume plots for four independent systems are within 2% deviation with respect to each other.

The measured average radius of gyration as function of temperature is not a smooth curve, therefore we can only consider the behavior of the chain length of 40 monomers.

The 〈R g 〉 value is bigger after the transition temperature compared to before the transition. (a) ( We further perform the cooling simulations for our big systems of different chain lengths from the temperature 350K to 250K at the cooling rate of 2K/ns using the Type-2B combined potential obtained from 40 monomers system. The volume per monomer as function of the imposed temperature is shown in Fig. 4.16. It can be noticed that the initial volume per monomer is higher for a short chain system due to the presence of more number of end monomers. The difference between big (2500) and small (20) system of chain length 40 monomers is that the drop in volume per monomer curve is much sharper compared to the small system. However, the transition temperature remains almost the same for both big and small systems. Unstrained NR crystallizes at a temperature between 293K and 223K [21]. Interestingly, the phase transition temperature in our model falls within this temperature window with the transition temperature increasing as function of the chain lengths. Note that for the big systems with 100, 000 beads, the DPD simulations are performed using the combined CG potentials. Bond length and angle distributions of the chain length of 400 beads system during cooling are shown in Fig. 4.17 where the distributions are measured in temperature windows of ∆ = 25K. Upon cooling the preference for larger angles increases as function of temperature. This indicates that segments of the chains prefer to form relatively straight sections because the angle potential has a minima at 180 • . In some polymeric systems that can crystallize upon cooling, for instance PVA, form distinguishable lamellae that are surrounded by amorphous regions, i.e., a semi-crystalline phase. In such a scenario the system undergoing cooling does not allow for the bond lengths to change much, even if the monomers are trapped within lamellae or amorphous regions or in between the ordered and the disordered regions. However, in our model the bond distance distribution becomes wider at lower temperatures and this observation is counter intuitive. To understand the behavior of our system during cooling, the radial distribution function or pair-correlation function g(r) is computed which provides information on the local positional order. In Fig. 4.18a the g(r) of all-monomers is shown and in Fig. 4.18b the g(r) of the inter chain monomers is presented for a long chain system at different temperatures. The first peak in the all-monomers g(r) corresponds to immediate bonded neighbors and the area under this peak reduces as function of a temperature which is in accordance with the observation in the bond length distribution. It is interesting to notice that the bond lengths, in the tail of the distribution are equal to the pair-distance of about 10 to 20 %. Consequently, the bonded particles of a monomer and the non-bonded particles will roughly have the same distance.

At high temperatures the curves show the expected liquid-like structure with regularly spaced broad peaks corresponding to successive neighbor shells. At low temperature (T = 250K) these peaks become more pronounced indicating some type of local order. To examine the structure in more detail we compute a perpendicular g ⊥ (r) and parallel g (r) correlation functions as mentioned in Chapter 3.3.3. The g ⊥ (r) is shown in Fig. 4.18c for inter chain monomers at different temperatures, whenever a pair of non-bonded monomers are within 6 Å. The peaks present at the low temperature suggest a two dimensional liquid like behavior of the non-bonded monomers.

In Fig. 4.18d, the g (r) correlation function is presented for the long chain system, whenever a pair of particles within 6 Å. We can observe from the curves at lower temperatures that the location of minima correspond to the integer multiples of the bond lengths, which indicates a preference for the alignment of neighboring chains. Based on these results we can say that at low temperature (250K) the aligned neighboring chains prefer to be offset with respect to each other.

The radius of gyration shown in the Fig. 4.19a is the average measurement per configuration as function of the imposed temperature. In the melt the 〈R g 〉 ranges between 15 and 50 Å for the chain lengths of 40 to 400 beads. Upon cooling the system undergoes a phase transition in the temperature range 270-285K and the size of polymers increases, i.e, the average R g changes to values between 30 and 80 Å.

The local bond order is computed using the largest eigenvalue S of the nematic order parameter 3.4. The same procedure is used to measure the orientational order as mentioned in Chapter 3. Note, in Chapter 3 crystallinity is defined as the fraction number of bonds with S ≥ 0.8. However, in this chapter we quantify the percent of locally ordered region by measuring the fraction number of bonds with S ≥ 0.6 as shown in Fig. 4.19b for different chain lengths. Note the choice S ≥ 0.6 is representative for increasing local order in a qualitative fashion. We can observe that the fraction of the locally ordered bonds increase with reducing temperature, which indicates the formation of some type of ordered region(s).

Snapshots of the melt and the semi-ordered configuration for a system of chain length 400 beads are shown in Fig. 4.20. The color coding used in this image is based on the

Anisotropy: CG system

It should be pointed out that the components of the pressure tensor become anisotropic below (nearly 10K) the transition temperature in our model. An example for this behavior as function of the imposed temperature for the long chain system is shown in Fig. 4.21a. The anisotropy of the pressure tensor arises because of the formation of bundle-like structure percolating across the simulation box during cooling. Whereas we use an isotropic pressure setting, for which, the diagonal components of the barostat are not independent of each other.

To avoid the anisotropic pressure tensor, we applied the anisotropic pressure setting in the Berendsen barostat where the 3 diagonal pressure components are not coupled. The simulation box adjusts to accommodate the imposed pressure in each direction independently. Using this method, we notice that the diagonal components of the pressure tensor remains nearly unchanged even at low temperatures which is well below the transition temperature of nearly 270K. However, the non-diagonal pressure components deviate from its average value of zero below the transition temperature as shown in Fig. 4.21b.

To overcome the deviation of the non-diagonal pressure tensor components, we used the triclinic pressure control using a Nosé-Hoover barostat. We switched from using the Berendsen barostat to the Nosé-Hoover barostat, as the standard module for the Berendsen barostat available in LAMMPS is does not de-couple the non-diagonal pressure elements Pxy, Pxz, and Pyz. Upon using the triclinic setting, all the components of the pressure tensor fluctuate around an average value as shown in Fig. 4.21c. The constant cooling simulation using the triclinic setting is performed from 280K to 230K at the same rate as mentioned before. The initial configuration at 280K is selected from the anisotropic pressure control simulation where the average pressure tensor components are at the desired value. Using the triclinic pressure control option we can cool the system to even lower temperatures than what is shown here. However, the shape of the simulation box becomes increasingly skewed.

The volume per monomer as function of the imposed temperature obtained using different options to control the pressure tensor elements are shown in Fig. 4.22a for the long chain system. This particular measurement along with the bond length and angle distributions are not affected by the pressure control methods. However, the average radius of gyration as function of the imposed temperature slight increases at low temperatures using the different pressure controlling settings in the increasing order starting from the isotropic, anisotropic and triclinic as shown in Fig. 4.22b.

Uni-axial Strain

As mentioned in the introduction 4.1 there are several experimental studies [22,28,147,[153][154][155]175] related to the strain-induced crystallization in both vulcanized, and unvulcanized forms of natural and synthetic rubbers. However, to our knowledge a molecular scale study of strain-induced/assisted crystallization on cis-PI has not been performed yet.

Having characterized the structure of the systems during a constant rate of cooling for different chain lengths, we use the same type of measurements to examine the effect of strain on those structures. Due to the finite size effect of our system as mentioned in 4.4.2 we consider the high temperature configurations. To this end, we perform an uni-axial deformation to the configurations at two temperatures, 350K and 300K for different chain lengths. We will focus on true strain rate of 10 -7 /fs, during which the lateral dimensions are independently subjected to atmospheric pressure and the temperature remains constant. The configurations are subjected to 50ns of NPT simulation before applying uni-axial strain. The results presented in this section are based on single realizations, but we have confirmed the qualitative behavior of the uni-axial straining by deforming the configuration in a different direction.

In Fig. 4.24, the stress-strain plots for systems with different polymer length are presented at two temperatures. The average stress is measured during the straining simulations over short time-period and the data have been fitted with a continuous piecewise function to guide the eye. The initial configurations are in the melt with weak local ordering. During initial straining the response of the stress is weak, followed by a faster increase of the stress beyond a certain strain. The easy identification of the strain beyond which faster increase of stress occurs depend upon the polymer length, the temperature and the strain rate.

If we focus on the system of short chain lengths of 40 and 80 monomers, the stress as function of imposed strain has a weak linear response for both the temperatures 350K and 300K. In the case of long chain systems, we can identify two types of stress behavior. Initially there is a weak linear response of the stress followed by a faster increase of stress due to the building up of stress. At temperature 300K the strain at which stress increases faster can be easily identified by visual observation compared to the high temperature of 350K.

It is important to point out that in the experimental studies [28], the vulcanized synthetic rubber exhibits the onset strain at around 4. However, in the case of unvulcanized IR such a behavior is detected at low temperatures of ≤ 273K, but the onset strain is reported in the same range. It appears that our CG model indicates the change in stress response at the same strain as observed in the experiments. Note that, for the chain length of 400 monomers, we have also performed the uni-axial deformation simulations at a slower strain rate of 10 -8 /fs and at a faster true strain rate of 10 -6 /fs at both temperatures. The stress-strain plot for the long chain system at different true strain rates are shown in Fig. 4.25. We can observe that for the slow strain rate the stress response to the imposed strain is similar to that for the short chain systems at 350K. During the slow strain rate the system has enough time to disentangle resulting in a weak linear stress response. However, during the faster strain rates the entanglements become more tight instead of disentangling, this results in building up the stress which can be identified as the domain followed after the initial linear stress behavior. During the fast strain rate the system has less time to adjust with changing strain values resulting in faster increase of stress at both the temperatures compared to slow true strain rates. At low temperature 300K the system has more local orientational order compared to high temperature therefore resulting in higher stress upon uni-axial deformation.

To quantify the entanglements present in the system we perform a primitive path analysis. The ends of the polymer chains are fixed. Next, each polymer chain is replaced by a series of particles that form impenetrable small straight segments. In other words, each bond that connects two consecutive monomers of the chains are replaced by several particles where the distance between the newly introduced particles prevents any possible bond crossings. Further, we perform a simulation which increases the tension between each particle of the same chain iteratively. To perform this simulation, we have changed the particle interactions which constitutes only the harmonic bond and non-bonded potentials. At the end of this procedure each multiple disconnected path becomes the shortest possible path. Since, the polymer-chain crossability is restricted, the number of locations where different polymers tries to cross each other is identified as an entanglement point. Note, we do not account for self entanglements therefore, due to this we underestimate the number of entanglements in the system. Another limitation of this procedure is that it can create entanglement at the ends of the chains because the end monomers are fixed therefore, the number of entanglements can be overestimated. A snapshot of a configuration with the primitive path of all polymers with chain length 400 monomers and an individual chain along with its corresponding primitive path are shown in Fig. To compute the entanglement length and the Kuhn length of the primitive path chain, we use the procedure described in Ref. [165]. After the relaxation of the configuration we can measure the chain length along the primitive path of a chain by summing up the bond lengths of the beads (L pp i ). Since the primitive path is supposed to obey Gaussian chain statistics, the Kuhn length is given as

a pp i = R 2 ee,i L pp i , ( 4.2)
where R ee,i is the end-to-end distance and L pp i is the chain length along its primitive path of the i-th chain. The average value of the entanglement length is computed using

〈N e 〉 = (N mono -1) 〈R 2 ee 〉/〈L pp 〉 2 , (4.3a) 〈R 2 ee 〉 = 1 M M i=1 R 2 ee,i , (4.3b
)

〈L pp 〉 = 1 M M i=1 L pp i , (4.3c)
where N mono is the number of monomers per chain and M is the number of chains in the system.

The results of the Kuhn length and the entanglement length for different chain lengths and temperatures are presented in Fig. 4.27. The value denoted as experiment in the plots 4.27b is obtained by a simple linear interpolation from the entanglement masses mentioned in Ref. [173] for two temperatures 298K and 413K. After interpolating the obtained entanglement mass is divided by the mass of a monomer (68.114 g/mol) to get the entanglement length in terms of number of monomers. The values obtained in our model for the entanglement length is in good agreement with the literature Ref. [162,165]. We also compute the entanglement length during straining for different chain lengths at 350K and 300K as shown in Fig. 4.28. We can see that the entanglement length increases as function of strain for long polymers (≥ 80 monomers per chain), indicating disentanglement occurring due to uni-axial straining. Note that for the chain length of 40 monomers, the average entanglement length increases as function of strain similar to the other chain lengths, but this increase in 〈N e 〉 is small compared to other chain lengths.

In the case of long polymers (≥ 400 monomers) the system has less entanglements at 350K compared to 300K. However, for the short chain systems (≤ 160 monomers) the 〈N e 〉 remains the same irrespective of temperature.

To investigate the structural measurement at the monomer scale, the bond length and angle distributions are measured during straining at the constant rate of 10 -7 /fs in 10 strain windows. The results of such measurements at a constant temperature 300K of the system with chain length of 400 monomers are shown in Fig. 4.29. Although the shift is small, the peak of the bond distribution shifts to higher lengths continuously as function of strain.

In the case of angle distributions, the population at larger angle increases with strain. This behavior becomes noticeable beyond the strain value 3. Similar to the case of cooling we can investigate the change in the structure using interchain monomer-monomer pair-correlation functions. For the long chain system at temperature 300K the results obtained are shown in Fig. 4.30. The same procedure is followed as mentioned in Chapter 3 (Section 3.4), where the strain direction forms as a reference for the perpendicular and parallel pair-correlations. The same observations that are made during cooling with respect to the local structure can be interpreted for straining systems as well, where the beads of the neighboring chains align and prefer to be shifted with respect to each other. The ordered regions during uni-axial deformation will align with the straining direction unlike the cooling where the ordered domains have random orientations. In the radial distributions Fig. 4.30c we can notice that for strain value of 3 the intermediate structure between melt-like at the low strain and the semi-ordered behavior at high strains (> 3) occurs in the regime where the stress increases sharply as function of the imposed strain in Fig. 4.24. Note that the RDF's are normalized to reach 1 at large distances. The fact that this does not appear to be the case is because the monomers of the same chain are not considered. To quantify the structure at monomer scale, we measure the local orientational order using the largest eigenvalue (S) of the nematic tensor Eq. 3.4. The results for the local ordering, i.e., the fraction number of bonds which are identified with S ≥ 0.6 are shown in Fig. 4.31 for different chains under constant true strain rate 10 -7 /fs for two temperatures. To quantify the amount of local ordering in the system we use the same definition during quenching and straining. It can be noticed that the amount of local ordering in our systems at low temperature 300K is nearly double compared to high temperature. We recall that the shortest chain system at low temperature (250K) has nearly 70 % of the bonds that have local orientational order S ≥ 0.6, but under staining they remain in a melt-like state. This can highlight the role of entanglements in the strain-assisted ordering. The short polymers have relatively few entanglements and upon staining these short polymers disentangle forming nearly straight chains that are oriented along the straining direction. Whereas, in the case of long polymers, the disentanglement time would be much longer compared to the simulation duration considered here. Therefore, the amount of local ordering is higher under straining at constant rate 10 -7 /fs when compared to the quenched system for long polymers. It is evident in the systems considered that the amount of locally ordered regions depend on temperature and the chain length. The effect of strain rates is illustrated by measuring the distribution of the local orientational order for three different true strain rates of the long chains of 400 monomers at 300K as shown in Fig. 4.32. During straining the peak of the distribution shifts to larger value for all the strain rates. In the case of fast true strain rate, the local orientational order distribution becomes double peaked at larger strain (≥ 4). However, in these strain windows the S distribution remains uni-modal that has longer tail for the slow strain rate (10 -8 /fs). In the case of intermediate strain rate of 10 -7 /fs, the S distribution exhibits a shoulder at strain ≥ 5. This indicates that the fast strain rate induces more local order in the system as compared to that of slow strain rates. The double peak in the local orientational order distribution typically represents the existence of different types of local ordering. However, to characterize the types of ordering induced in the system due to fast strain rate, more detailed analysis is required.

Another type of local order is the second Legendre polynomial of the cosine of the angle between the unit bond direction b and the unit strain direction ε as mentioned in Eq. 3.5. The results of this measurement are presented in Fig. 4.33 for different chain lengths at two temperatures. This measurement only indicates that the average bond director increases along with strain which is expected, as we are deforming the simulation box in one direction. In the case of short chains, the P2 measurement appears to reach a plateau at large strains. This indicates that the short polymers have a liquid-like behavior. In the case of long polymers (≥ 400 monomers) this P2 measurement appears to increase continuously along with strain. Whether this measurement continues to increase upon further straining or reaches a plateau can not be determined from this study.

Based on these structural measurements we can identify reorientation of the bonds and in turn polymers occur in the first regime, i.e., at short strain where stress-strain behavior has a weak linear relationship. Next, the stretched polymers or chain segments get aligned with respect to each other. Further, we examine the structural changes taking place during the uni-axial deformation at the polymer scale by measuring the radius of gyration tensor using the Eq. 3.6. Since the simulation box is deformed along a particular direction, we use the components of the R g tensor along the parallel and the perpendicular straining directions as function of the imposed strain to characterize the shape/size of the polymer as shown in Fig. 4.34. As expected, the strain on the system stretches along the straining direction and compresses the polymers in the perpendicular straining direction. To quantify the structure at the polymeric scale, the end-to-end distance is shown for different chain lengths at two temperatures in Fig. 4.35 as function of the imposed strain. way to improve the CG model with (Type-2 variants) is to include more than one type of CG bead per chain that corresponds to the following:

• Middle bead: CH CH 2 CH 2 C(CH 3 )

• End bead-1: CH 3 C CH 3

• End bead-2: CH CH 3

By using such a mapping scheme the density related measurements can be better reproduced in the CG model compared to the simple mapping scheme we have used.

Another way to improve the CG model is to use a torsion potential which can help in better reproducing structure at polymer scale like the radius of gyration and end-to-end distance.

One more way for improve our model is to use piece-wise polynomials instead of tabulated potentials. In the range where the occurrence is limited in the tabulated potentials, there will be considerable uncertainty which might result in undesirable forces.

In this study we have used the settings based on simpler models [80,110,111] and assume that these settings are still valid. This can be another alternative approach to further improve our CG model which is fine tuning the cut off values (conservative and dissipative forces).

We have assumed that the CG potentials are temperature independent in this work. But this is not true so to obtain a realistic CG model we have to generate a temperature dependent coarse-grained potential. To do this we need to perform atomic simulations at different temperatures. As we go to low temperatures, the atomic simulations will take considerably long time to equilibrate which could well be beyond our computational resources. To solve such an issue, back-mapping can be used, where the CG systems can be used to obtain the atomic system as mentioned in Ref. [176].

Having mentioned the possible ways to improve the simple (Type-2B) CG model, we summarize the findings as follows. For short chain lengths of 10 and 20 monomers in the Type-2B CG model, the effect of end monomer is predominant compared to other CG models (Type-1A and 2A). The long chain system of 80 monomers suffers from slow dynamics which affects the CG potentials and hence CG models of 30 and 40 monomers are equally good compared to other chain lengths. We mention again that for the chain length of 40 monomers four independent realizations were each subjected to STM process and these four CG potentials were combined to get one coarse-grained potential. Using the angle and bond length distributions of the combined potential of 40 monomers we generate bigger and longer chain systems.

These systems were subjected to continuous cooling at various constant rates where we focus on 2K/ns. We can see that the system undergoes a phase transition in the temperature range 270K to 285K. It is interesting that the unstrained natural rubber crystallizes approximately in the same temperature window (293K to 223K) [21]. The structural quantities at various scales ranging from macroscopic measurements like temperature-volume plot to monomer scale like pair-correlation functions are measured for different chain lengths while cooling from 350K to 250K. It is evident from various structural measurements at different scales that the transition from melt to semi-ordered structure occurs using this CG model (Type-2B). From the perpendicular and parallel pair-correlation functions we can see that the alignments of neighboring chain occurs at low temperature (250K) but not the formation of distinguishable crystallite (lamellae). The visual inspection also suggests the same. During continuous cooling we encountered a pressure anisotropy in our system where the pressure tensor components upon cooling become anisotropic. To overcome this issue, we used different barostat setting like de-coupling the pressure components using anisotropic and triclinic options. Using the triclinic barostat setting the pressure tensor components remain isotropic even at low temperature below the transition temperature. The use of different barostat options does not affect some measurements like the transition temperature, but at the polymer scale this effect can be seen as the average radius of gyration increases from isotropic to anisotropic and then to triclinic barostat setting at low temperature.

The strain induced or assisted crystallization is observed for long polymer systems. During the initial straining, polymers are quite flexible which allows for reorientation of the bonds. This is followed by the local alignment of the chains along the straining direction. We have performed the uni-axial straining study for different chain lengths at two temperatures. We have shown as a representative, the behavior of 400 monomer system at different strain rates at both temperatures of 300K and 350K. We have measured the entanglement length for the unstrained system at both temperatures. The measured entanglement length computed in our model is in good agreement with the literature finding for other CG models and experiments [162,165]. Upon performing a similar analysis for the strained configuration we find the entanglement length increase as function of strain for different chain lengths. From the local structure analysis we find that the preference for larger angles increases as function of strain. From the pair-correlation functions, we observe the local ordering is similar in both cooling and straining simulations, the alignment of the neighboring chains or chain segments occur. We find local ordering similar to that of bundle-like structure but not the formation of lamellae in our CG model.

Chapter 5

Network Systems

Introduction

Vulcanized NR has been studied widely because of its superior mechanical properties compared to those of synthetic rubbers like isoprene rubber, polybutadiene rubber, and styrene butadiene rubber. The strain induced crystallization of rubber is considered as a reinforcing mechanism to prevent premature rupture and slow down crack propagation [177].

Flory [17] proposed a theory to describe strain induced crystallization for a system with high molecular weight like vulcanized rubber. He described the stress-strain relationship to be weak at small strain which is followed by a steep increase of stress as the strain increases. He associated the steep increase of the stress with the decrease of accessible chain configurations when the system is strained, as chains between cross-links become close to fully extended.

Non-Gaussian chain network models with the classic theory of rubber elasticity have been used to study the stress response upon straining of sulfur vulcanized NR [178,179]. Arruda and Boyce [164] successfully simulated the stress-strain behavior of the experimental data presented by Treloar [180] under different loading scenarios including uni-axial deformation. The success of reproducing the stress-strain relation of sulfur vulcanized NR in these studies suggest that SIC has no contribution to the stress [28].

In the above mentioned works of Flory [17] and Treloar [180] cross-links were treated as pivot points in a network. Following these seminal works, several experimental studies on strain induced crystallization in rubbers have been reported [12,22,[181][182][183]. In particular, real-time X-ray diffraction experiments have been used to reveal the ability of different materials, including rubbers, to crystallize during straining. Many studies have used a synchroton X-ray source associated with a high speed stretching machine to estimate a characteristic time of SIC around 50 to 100 milliseconds by means of wide angle X-ray diffraction [177].

Toki et al. [28] performed uni-axial tension on both vulcanized and unvulcanized samples of NR and IR at various temperatures. They reported that at high temperature vulcanized IR shows SIC which is not observed in unvulcanized IR. On the other hand, SIC is observed in both vulcanized and unvulcanized NR at high temperatures. The presence of pseudo end-linked network which is reported to be constituted of non-rubber components in NR has been raised as an explanation for this observed difference [153].

Sainumsai et al. [177] also performed tensile tests of both unvulcanized and vulcan-ized NR and IR for different cross-link densities. The stress-strain curves obtained for unvulcanized-NR, vulcanized-NR and vulcanized-IR are shown in Fig. 5.1. They perform tensile straining experiment at 297K with the initial strain rate of 0.5s -1 .

Figure 5.1: The stress-strain relations for rubbers. The uni-axial straining is performed at 297K. Image taken from Ref. [177].

Sainumsai et al. report that during straining the stress increases linearly until a certain strain. Beyond this strain, the behavior of stress becomes different from the first regime which is referred as stress-upturn and the corresponding strain is known as the onset strain. Their results indicate that the onset strain signals strain induced crystallization of vulcanized NR and is independent of the cross-link density. However, in the case of vulcanized IR the onset strain appears to be dependent on the cross-link density.

Different vulcanizing agents and curing processes are used for obtaining vulcanized rubbers. In the experimental studies mentioned above [28,177] sulfur is used as a vulcanization agent and different ratios of cross-link structures (monosulfide, disulfide and polysulfide) can be obtained. More details on the preparation of vulcanized samples can be found in Ref. [177].

The network systems have also been investigated by various simulation methods to probe the different length and time scale behaviors associated to these systems. Grest and Kremer [184] used MD simulations to generate polymer networks in the early 1990s where cross-linking was done on well equilibrated linear polymers represented by a bead-spring model with LJ interactions. They investigated the statistical properties of random cross-linking in polymer melts. The cross-links were placed at random in the system and a bond was created between the cross-link particle and the closest monomer. This monomer represents an "activated" particle in the system. The second monomer that connects to the cross-link particle was chosen based on a distance criteria.

The first fully atomistic MD simulation of a polymer network with less than 200 atoms was studied by Hamerton et al. [185]. They predicted elastic modules and T g for trifunctional polycyanurate which were in reasonable agreement with experimental data. In the late 1990s Doherty et al. [186] performed progressive polymerization reactions using MD simulations. The cross-linking process involved, reactions of monomers with each other based on distance criteria. The formation of chemical bonds was based on the distance between pre-defined reactive sites on the monomers and growing of chains. Their crosslinking procedure involved packing and equilibrating monomers, repetitive cross-linking until a particular conversion is achieved and the relaxation of final structures based on energy minimization. Several other computational procedures for constructing network systems have been explored for various polymer systems. More details on network formation can be found in a review article [187].

Edwards et al. [188] developed a slip link model to study the dynamics of cross-linked polymer networks. In this model, the motion of polymer segments is considered to be confined in a tube-like section and the entanglements are treated as discrete but mobile constraints. More recently, Munoz et al. [189] 1 used elastic network model where the cross-links were treated as mobile constraints to study the stress-strain response in a cis-PI based system.

At the molecular scale, Li et al. [162] generated cross-linked systems based on a cis-PI coarse-grained model. They used a CGMD simulation with a distance constraint between cross-linked PI beads and corresponding sulfur beads. The sulfur beads in their model represents two sulfur atoms. They introduced the cross-links randomly and if the distance between a cross-link bead and the cis-PI monomer is within a certain distance bonds were formed between the cis-PI and the cross-link beads. They did not allow the cross-linking within the same chain. Using the network model, they performed a primitive path analysis on different cross-linked systems.

Here we study the stress-strain behavior of a network system using a cis-PI inspired model. To our knowledge a molecular scale study of strain-induced/assisted crystallization on cis-PI network system has not been performed yet.

Methodology

In the process of obtaining the CG potentials as mentioned in Chapter 4 (Section 4.2) a full atomistic model for polyisoprene is simulated using the AMBER force field. Four independent atomic models of chain length 40 monomers having 20 polymers are each subjected to 150ns of equilibration run at 350K. Next, MD trajectories are obtained where a configuration is stored every 50fs. The atomic trajectories are used with the Type-2B coarse-grained mapping to get the reference trajectories. In the Type-2B CG mapping scheme, the center of the CG bead is located at the center of carbon-carbon single bond as depicted in Fig. 4.1c. In this CG mapping choice the end atoms are neglected which effectively means there is one monomer less in the coarse-grained representation of a polymer. The reference trajectory is subjected to the STM process by selecting the following interaction styles: harmonic bond, tabulated angle and pair-wise interactions. This is done for all the four systems independently and then one combined potential is obtained. The interactions of the combined potentials are shown in Figs. [START_REF] Biron | Material Selection for Thermoplastic Parts[END_REF].5(e, f) and 4.6c that correspond to harmonic bond, tabulated angle, and tabulated pair-wise interactions, respectively. These CG potentials are given in Appendix A.2. Note that we have used the same interactions to perform cooling and straining simulations for big systems with 100, 000 beads in Chapter 4.

In this chapter we perform uni-axial deformation of network configurations using DPD simulations. We use the big systems with 100, 000 beads of different chain lengths 40, 80, 160, 400 and 800 monomers to generate polymer networks. The procedure of generating these big systems was discussed in Chapter 4 (Section 4.4.1). The equilibrated configurations are used to generate network systems with different cross-link density. The cross-link density is computed with respect to the total number of monomers (cis-PI beads) in the considered system. For instance 1% cross-link density corresponds to 1, 000 cross-link beads for a system with total of 100, 000 monomers.

To generate the network system, we first determine every pair of monomer in different polymers or separated by at least 4 bonds in the same chain and a distance between 2 and 6 Å. Then a cross-link bead is introduced half way, i.e., at the mid point between every monomer pair selected. Later, from this complete set of possible cross-links the desired number is selected randomly such that a monomer has at most one cross-link and the crosslink bead is separated by more than 2 Å from every other cis-PI or cross-link bead.

In this model we can identify three types of beads:

• Type-1 is the cis-PI monomer bonded with two other cis-PI monomers of the same chain

• Type-2 is the cis-PI bead that has three bonds where two bonds are connected to cis-PI monomers of the same chain and one bond connected to the cross-link bead

• Type-3 is the cross-link bead that represents a disulfide monomer which has two bonds connected to two cis-PI monomers Note that we can also identify another type of bead, end-monomer which has only one bond, since all the interactions are same for a middle and end monomers, we classify the end beads as Type-1 bead.

Here we take a simpler approach to generate the CG network systems and apply the combined potentials obtained in Chapter 4 for a cis-PI system. The interactions for crosslink beads are the same as those of cis-PI monomers. The angle interaction for the cross-link beads is turned off. Only the bond interaction, i.e., the harmonic bond, and the non-bonded interaction of the combined potentials are applied to the cross-link beads.

The generated network systems are subjected to a NPT simulation at a constant temperature of 350K and constant pressure of 1 atmosphere. We have used the Berendsen barostat with anisotropic setting and the simulation is performed until the density of the system converges. The measured densities for different chain lengths with different cross-link densities are shown in Fig. 5.2. The mass density of the systems increases with the increase in cross-link densities for all chain lengths similar to uncross-linked system as the effect of end particles is higher for short chains. Beyond the chain length of 160 monomers the mass densities of the systems remains nearly the same irrespective of the cross-link density. Later, these systems are subjected to uni-axial deformation at a constant true strain rate during which the lateral dimensions are independently controlled by a barostat set to atmospheric pressure and the temperature remains constant. Here, we focus on the true strain rate of 10 -7 /fs with different cross-link densities at a constant temperature of 350K.

= 0 % crosslink = 1 % crosslink = 2 % crosslink = 5 % crosslink = 10 % crosslink = 20 % Expt MD: crosslink = 0 %

Uni-axial Strain

Here we perform a similar analysis for the network systems as done for straining the melt systems presented in Chapter 4 (Section 4.4.3). The bond length distributions for systems with different cross-link densities of chain length 400 monomers are computed at different strain windows as shown in Fig. 5.3. The bond lengths do not change much during uni-axial deformation for non-network system, as expected, because it is mainly the angle potential that is responsible for chain elongation along the straining direction. We can observe this from the angle distribution 4.29b in Chapter 4. Note that the angle distributions for the network systems are similar to that of the non-network system.

The bond distributions corresponding to larger strain windows of the network systems show bigger tails which accounts for longer bond lengths. This is a consequence of the harmonic bond potential, which allows for over-stretching of the bonds. However, in this study we make a simple choice by selecting a hard cutoff for the bond length and all the measurements beyond this cutoff are not considered. We recall that at the atomic scale, beads are composed of four carbon atoms in the chain backbone. The fully stretched backbone would roughly correspond to 5.5 Å. Therefore, we do not consider the measurements when more than 1% of the total fraction of the bond distances are larger than this cutoff distance. This means for higher cross-link density systems we effectively have only a small strain range where our model has a physical behavior.

The stress-strain plots for different chain lengths with 1% cross-link density are presented in Fig. 5.4a. In the case of short polymer systems with a chain length of 40 monomers, the stress appears to reach a plateau similar to that of non-network system. Since these systems correspond to low cross-link density where on average we expect to find less than one cross-link per chain, it behaves very similar to a non-network system. In the case of longer chains (≥ 80 monomers per chain), the stress response remains similar up to a certain strain with a linear behavior. Beyond this strain, the effect of cross-links starts to dominate and as a result the stress increases more rapidly. The divergent behavior of stress is mainly because of harmonic bond potential. This phenomena is predominant for long chains (≥ 400 monomers) compared to short polymers (80 and 160 monomers), as these systems on average will have 8 to 16 cross-links per chain for 400 and 800 monomer chain length systems, respectively. .4b depicts the stress-strain plots for the 400 monomer chain system with varying cross-link densities. It can be noticed that for cross-link density ≥ 2%, the strain range is truncated based on the bond length cutoff criteria. It can be observed that the higher the cross-link density smaller is the strain range where the linear stress response occurs. This initial domain is followed by the divergence of the stress, i.e., the non-linear stress-strain behavior. The divergence of stress after reaching a particular strain is observed in many experimental works for vulcanized IR [177,182,190]. When the system has more crosslinks, a path between a monomer and its image over the periodic boundary condition will exist, i.e., a percolating network is formed. For higher cross-link density systems, the system realizes it is a percolating network at smaller strain this results in the shorter initial regime where the stress response is linear.

To analyze the local structure at monomer scale we compute the local orientational order. We use the nematic order tensor Q αβ expressed in Eq. 3.4 and calculate the largest eigenvalue S. The local orientational order S is computed by defining a local spherical region for every bond in the system, where the position of the bond is taken to be the mid point between the corresponding monomer pair and the sphere contains on an average of 14 particles similar to the PVA and the cis-PI systems in Chapters 3 and 4, respectively.

φ can be used to quantify the amount of local ordering by using the same definition mentioned in Chapter 4, where the fraction of bonds (S ≥ 0.6) is considered. This measurement as function of the imposed strain for different chains with 1% cross-link density is shown in Fig. 5.5a and in Fig. 5.5b the same quantity is presented for the long polymers of 400 monomers per chain with varying cross-link densities. As expected the amount of local ordering increase with increasing strain for all the chain lengths irrespectively of the cross-link densities. Note that the selection for S, i.e., the cutoff value for φ is based on the distributions of the local orientational order measured in different strain windows like those shown in Chapter 4 (Fig. 4.32). We can notice that the 〈S〉 increases systematically with strain. Therefore, the choice made for cis-PI model in this Chapter and in Chapter 4 for the fraction bonds with S ≥ 0.6 is representative for increasing local order in a qualitative fashion. We find a similar behavior with φ(S ≥ 0.5), φ(S ≥ 0.7) or any other cutoff. Overall, the amount of local ordering will change along with the "onset" strain value, i.e., the strain value at which the fraction of the bonds starts to increase by leaving the initial region, but the qualitative picture still holds true. Another type of local order measurement is the second Legendre polynomial Eq. 3.5. This measurement provides the information on the average bond direction along the strain direction as shown in Fig. 5.6 as function of the imposed strain for different chain lengths and different cross-link densities. This highlights the reorientation of the bonds at short strain followed by a gradual orientation of the polymers along the straining direction. In the case of network systems, the polymers cannot continue to align due to the presence of permanent bonds between different chains which restricts the alignment.

To examine the structural changes at the polymer scale we present the results of the radius of gyration components and the end-to-end distance as function of strain as shown in Figs. 5.7 and 5.8, respectively. We measure these quantities for every chain and then average it over all the polymers present in the configuration. The cross-link beads are not considered for any of the structural measurements and therefore, an individual polymer is a composed of only cis-PI beads. The radius of gyration along (R ,g , solid lines) and perpendicular (R ⊥,g , dotted lines) to the strain direction as function of the imposed strain at a constant temperature 350K for (a) different chain lengths with 1% cross-link density and (b) chain length of 400 monomer system with different cross-link densities under the true strain rate of 10 -7 /fs.

As is to be expected, the strain on the system stretches the polymers along the straining direction and compresses the polymers in the perpendicular strain direction. In the case of long polymers (≥ 400 monomers) the average R g components reach approximately the same value at large strain for both cross-linked and uncross-linked system. Similar observations can be made for the average end-to-end distance for the different types of system.

Summary and discussion

We have generated the network systems of different chain lengths with varying cross-link densities at 350K. The cross-links have the same bond and pair-wise interactions as the cis-PI monomers. The angle interaction for cross-link beads is not included. We have taken a simple approach to generate the CG network. However, a better way to simulate a coarsegrained network is to generate an atomic configuration with cross-link atoms (like sulfur atoms). After equilibration of the atomic network, the MD trajectory is obtained where a configuration is stored at regular interval. Later, the MD trajectory should be subjected to geometric coarse-grain mapping to get the corresponding reference trajectory. After which the appropriate interaction styles must be selected and then subject to the STM process to get the CG potentials. The choice of bond interactions must be selected bearing in mind the problem associated with bond stretching that will arise in a CG network during straining. One solution could be to use finite extensible nonlinear elastic (FENE) potential and another way to perform a realistic cis-PI network simulation at the molecular scale is to allow for bond breaking in some fashion that will prevent high bond energies.

The mass densities of the system increase with increasing cross-link densities for all the chain lengths which is similar to the molecular length dependent behavior as exhibited for uncross-linked system. The network system is then subjected to a uni-axial deformation at different constant true strain rates where the results of (10 -7 /fs) true strain rate are presented. The bond length distribution measured at different strain windows shows that for high (≥ 2%) cross-link density systems, an extensive bond stretching will take place. This is the consequence of using harmonic bond interactions. To overcome this unphysical behavior, we put a hard limit on the bond length where we do not consider any measurements when the system has more than 1% of the total fraction of the bonds larger than a 5.5 Å cutoff length.

Overall the picture that emerges for long polymer systems at temperature 350K is that the reorientation of bonds and alignment of chains occur due to uni-axial straining. The initial response of the system is the reorientation of the bonds that occurs due to the stretching of the polymers along the straining direction. In the next regime the system realizes that it is a percolating network. The deformation of this network structure results in a non-linear (divergent) stress response. During this stage, the system becomes more locally ordered but the individual chains cannot continue to fully align due to the permanent connection between different chains. At the polymer scale the effect of cross-links is not much as the average end-to-end and the radius of gyration components reach a similar value at large strain. We make similar observations with respect to the stress response, local order and polymer scale measurements as function of strain for slow (10 -6 /fs) and fast (10 -8 /fs) true strain rates. The main difference being the 'onset strain' after which there is a faster increase stress and other measurements changes for different true strain rate simulations. Qualitatively the observations made for 10 -7 /fs true strain apply to other true strain rates as well.

Chapter 6

General conclusions and perspectives

The SIC of polymeric systems is a pretty old topic. The structural properties and their associated relationship with SIC is well investigated from a macroscopic point of view. Crystallization in polymers can be induced via temperature, i.e., cooling from melt, mechanical stretching or solvent evaporation. The crystallization of polymers affects mechanical, optical, thermal and chemical properties. The crystallization process is associated with the partial alignment of chains forming a semi-crystalline phase. Since more than two-thirds of the world's polymeric applications are covered by semi-crystalline polymers, this field of polymer research has obtained great interest. Although SIC is well investigated using experiments, there are certain aspects that are unclear. For instance, germination of crystallites and how do these crystalline domains change upon deformation at the scale of crystallites. The rising industrial constraints like recyclability, biodegradability and bio-sourcing puts forward new challenges to the synthesis of novel polymers exhibiting SIC. Such challenges require the ability to anticipate the structural and kinetic properties at the proper mesoscopic scale while retaining the most relevant aspects of the chemical structure at an atomic level in order to direct new molecular design.

It would be ideal to study any polymer system at high resolution, but crystallization is a multi-scale process with length-scales ranging from Ångströms to micrometers and timescales from femtoseconds to milliseconds. This has influenced the development of different coarse-grained methods.

We have focused on a bottom-up approach where the information from the atomistic scale is passed to the mesoscopic scale using the STM procedure for cis-polyisoprene. To obtain the CG model for cis-PI atomic systems, we perform MD simulations for different chain lengths at 350K using the AMBER forcefield. Next, we apply different CG mapping schemes and obtain the CG potentials using the STM process. The DPD simulations are performed using the obtained CG potentials and the results are compared with the reference MD. We identify a particular CG mapping scheme (Type-2B) and obtain the combined CG potentials based on four independent systems.

The CG model for PVA is obtained using the iterative Boltzmann method. The coarsegrained MD simulations are performed using the interactions obtained from literature [123,124].

We perform temperature-induced and strain-induced/assisted crystallization studies of both the CG models, PVA and cis-PI systems. The studies pertaining to temperature-induced crystallization is done by cooling the melt system continuously at a constant rate and SIC is done by means of uni-axial deformation. To characterize the structure we have developed various measurement tools to quantify the structure formed at different length scales.

Upon continuous cooling of the melt configurations at different constant rates for both the CG models, we find a phase transition to a semi-crystalline/ordered phase. In the case of CG-PVA systems we observe the formation of different crystalline domains below the transition temperature. However, in the case of CG cis-PI model we find the formation of bundle-like structure below the transition temperature, i.e., the co-existence of a continuous ordered structure accompanied by disorder regions but without any distinguishable lamella. To characterize the semi-crystalline domains formed in the CG-PVA system, we perform a cluster analysis based on the local orientational order.

Uni-axial deformation is performed at different constant true strain rates for both the CG models. In the case of CG-PVA model, SIC is observed in three stages. In the first stage, polymers are still flexible and mobile thereby allowing for the local reorientation of polymer bonds. In the second stage, rapid crystallization was seen and at the last stage, the local structures, such as hair-pins, were removed from the system obtaining almost fully elongated chains. In the case of CG cis-PI model, the strain induced or assisted crystallization is observed as a continuous process. During the initial straining, polymers are quite flexible which allows for reorientation of the bonds. This is followed by the local alignment of the chains along the straining direction. In the case of network systems we use the CG potentials inspired from the cis-PI system where the initial response of the system is the reorientation of the bonds due to the stretching of the polymers along the straining direction that is similar to the uncross-linked system. In the next regime the system realizes that it is a percolating network. The deformation of this network results in a non-linear (divergent) stress response. During this stage, the system becomes more locally ordered but the individual chains cannot continue to fully align due to the permanent connection between different chains.

It is observed that the angle potential is predominantly responsible for chain elongation or formation of straight segments in both the CG models.

The CG model for the cis-PI systems is a first step to obtain a realistic CG potential using the STM process. Here we apply the STM procedure for the first time on a full atomistic representation to obtain the CG potentials. Note that we have made a few assumptions like that our CG model is temperature independent and that the settings (cutoff distance for dissipative and conservative forces) based on simpler models are still applicable for our model. Despite these assumptions our CG model is capable of reproducing some measurements like the transition temperature, density and qualitatively reproduces the stress-strain behavior well with respect to the experimental findings. With these promising aspects it would be worth-while to invest more resources to develop more realistic CG potential of cis-PI system.

Another limitation is concerning the STM approach itself. Under the current setup, the matrix operations work reasonable well for a small system of about 800 monomers. However, if we were to back map a 100, 000 monomer CG system to get an equivalent atomic configuration then such an atomistic system will have 13 times more particles. The STM process on such a big system would take considerably more computational resources to obtain a CG potential. In future, it might be a good idea to improving this aspect of STM to get realistic CG potentials based on the big atomic systems.

Modelling strain-induced crystallization at the molecular scale

Strain-induced/assisted crystallization (SIC) of polymers plays an important role in a variety of industrial applications, such as tires, fibers, and containers. The semi-crystalline structure that is induced by straining is assumed to be the origin of the high tensile strength and resistance against the formation and propagation of cracks in rubber. The process of SIC has been well investigated experimentally and remains a topic of active research. However, some of its aspects, such as the germination of crystallites and how these crystalline domains change under deformation at the scale of crystallites, are still not fully understood. The aim of this thesis is to investigate SIC of polymers at the molecular scale by means of computer simulations.

We have used coarse-grained (CG) simulations for two different systems. The first is a CG model for polyvinyl alcohol that is well established in the literature and was developed with the aid of the iterative Boltzmann method. The second CG model is for cis-polyisoprene and was obtained here (Chapter 4) from atomic simulations in combination with a Bayesian optimization method in an approach known as Statistical Trajectory Matching.

We investigate two types of crystallization, temperature-induced by continuous cooling of a melt system at a constant rate and SIC by means of uni-axial deformation. Both CG models are able to reproduce some measurements, like the transition temperature and density, and qualitatively illustrate the stress-strain behavior in accordance to what is observed experimentally.

This work inspires and encourages the development of a more realistic CG model for cispolyisoprene using the same multi-scale simulation method. These CG potentials can be used to examine much longer time and length scales that are inaccessible by atomistic simulations and provide new insights in the different phenomena that play a role in the process of polymer crystallization.

Modélisation de la cristallisation sous tension à l'échelle moléculaire

La cristallisation sous tension (CST) des polymères joue un rôle important dans une variété d'applications industrielles, telles que les pneumatiques, les fibres et les récipients. La structure semi-cristalline induite est supposée être à l'origine de la grande résistance à la traction, à la formation et à la propagation de fissures dans le caoutchouc naturel. Le phénomène de la CST a été étudié expérimentalement et reste à l'heure actuelle un sujet de recherche actif. Cependant, certains de ses aspects, tels que la nucléation des cristallites et la façon dont ces domaines cristallins évoluent sous déformation à l'échelle des cristallites, ne sont pas encore totalement compris. L'objectif de cette thèse est d'étudier la CST des polymères à l'échelle moléculaire au moyen de simulations numériques.

Nous avons utilisé des simulations gros grains (GG) pour deux systèmes différents. Le premier est un modèle GG pour l'alcool polyvinylique qui est bien établi dans la littérature et a été développé à l'aide de la méthode itérative de Boltzmann. Le second modèle GG concerne le cispolyisoprène et a été obtenu ici (chapitre 4) à partir de simulations atomiques en combinaison avec une méthode d'optimisation bayésienne connue sous le nom de Statistical Trajectory Matching.

Nous avons étudié deux types de cristallisation, celle induite par la température avec un refroidissement et continu d'un système fondu à une vitesse constante et celle induite par la CST au moyen d'une déformation uni-axiale. Les deux modèles GG sont capables de reproduire certaines mesures, comme la température de transition vitreuse et la densité, et illustrent qualitativement le comportement contrainte-déformation en accord avec ce qui est observé expérimentalement.

Ce travail inspire et encourage le développement d'un modèle GG plus réaliste pour le cispolyisoprène en utilisant la même méthode de simulation multi-échelle. Ces potentiels GG peuvent être utilisés pour examiner des échelles de temps et de longueur beaucoup plus longues qui sont inaccessibles par les simulations atomiques et fournir de nouvelles perspectives sur les différents phénomènes qui jouent un rôle dans le processus de cristallisation des polymères.

Figure 1 . 2 :

 12 Figure 1.2: Schematic diagram of volume per monomer as function of temperature showing different cooling paths at a constant pressure.

Figure 1 . 3 :

 13 Figure 1.3: An atomic representation of short branched low-density polyethylene. Image taken from Ref. [5].

Figure 1 . 4 :

 14 Figure 1.4: A representation of vulcanized polyisoprene.

Figure 1 . 5 :

 15 Figure 1.5: A schematic representation of (a) an elastomer network along with (b) crosslinks and (c) entanglements.

Figure 1 . 6 :

 16 Figure 1.6: A schematic representation of modulus-temperature relation for a polymer.

Figure 1 . 7 :

 17 Figure 1.7: Illustration of different length-scales found in the crystalline structure of polymers, ranging from the atomic to the spherulite level. The different images correspond to (a) atoms of a monomer, (b) collection of several monomers representing a polymer, (c) different polymers with folds and stacked together forming a lamella and (d) several lamellae connected by amorphous sections to form a spherulite (Image from Ref. [26]).

Figure 1 . 9 :

 19 Figure 1.9: Schematic diagram of chain-folding models. (a) Adjacent re-entry model, (b) Random-switchboard model and (c) Intramolecular nucleation model.

Figure 1 . 10 :

 110 Figure 1.10: Schematic representations (2D) of growth of different regimes according HLtheory. Each square corresponds to the cross-section of a stem.

Figure 1 . 11 :

 111 Figure 1.11: Schematic representation of a possible structure during rough crystal growth and frustration; the minimum crystal thickness (L c,min ) at a particular crystallization temperature. Image taken from Ref. [25].

Figure 1 . 12 :

 112 Figure 1.12: Schematic representation of coarse-graining steps. Image taken from Ref. [39].

Figure 2 . 1 :

 21 Figure 2.1: Geometry of a model short molecule representing bond distance r i j , angle θ jih , and the torsional angle φ k jih .

M

  o l e c u l a r s i m u l a t i o n s D P D / M D C o n t i n u u m m e c h a n i c s Q u a n t u m m e c h a n i c o r m a t i o n C G p o t e n t i a l o b t a i n e d f r o m S T M / I B M F o r c e f i e l d

Figure 2 . 3 :

 23 Figure 2.3: Hierarchical multi-scale computational approaches to predict the viscoelastic properties of polymers. Various time and length scales associated with different computational methods. Figure reproduced from [94].

Figure 3 . 1 :

 31 Figure3.1: The bending potential for the CG-PVA model[123,125]

Figure 3 . 2 :

 32 Figure 3.2: Schematic representation of an artifact formed after generating the initial configuration where the three polymers are shown in different colors with each having two bonded beads. The six beads form an octahedral structure.

Figure 3 . 3 :

 33 Figure 3.3: The average volume per monomer at different quenching rates Ṫ for short polymers, 100 monomers (solid lines), and long chains, 1000 monomers (dashed lines). The latter curves have been shifted upwards by 0.00269/σ 3 , the dotted curve is an unshifted reference. (Error bars are indicated only for short chains but are very small.)

Figure 3 . 4 :

 34 Figure 3.4: The average volume per monomer for long chains of 1000 monomers at different quenching rates (dotted lines) and heating (solid line).

Figure 3 . 5 :

 35 Figure 3.5: A sequence of conformations of a long chain at the cooling-rate Ṫ = 10 -6 τ -1 for the temperatures as indicated.

Figure 3 . 6 :

 36 Figure 3.6: The average radius of gyration per polymer at different cooling rates for (a) short polymers of 100 monomers and (b) long polymers of 1000 monomers

64 Figure 3 . 8 :

 6438 Figure 3.8: Probability density distribution P R g of the radius of gyration for polymers of length, 100 and different temperature windows during the cooling rate of Ṫ = 10 -6 τ -1 . The two insets are typical examples of polymer conformations corresponding to a single and double hair-pin structure observed for R g ≈ 7σ and 5σ, respectively.

Figure 3 . 10 :

 310 Figure 3.10: The average length-weighted stem length as a function of temperature for different cooling rates for the chain lengths of 100 (solid) and 1000 (dotted).

Figure 3 .

 3 Figure 3.11: (a) The radial distribution function g (r), (b) the two-dimensional perpendicular distribution g ⊥ (r), and (c) the one-dimensional parallel distribution g (r) for the long polymers at different temperatures during the slowest cooling rate. The parallel and perpendicular directions are with respect to the local chain orientation.

Figure 3 . 12 :

 312 Figure3.12: The fraction crystallinity χ (S > 0.8) as a function of the imposed temperature at different cooling rates for polymers of 100 (solid) and 1000 (dotted) monomers, respectively.

Figure 3 .Figure 3 . 14 :

 3314 Figure 3.13: (a) The average number of crystallites per 10 5 monomers and (b) the average size of cluster during cooling for polymers of 100 (solid) and 1000 (dotted) monomers, respectively.

Figure 3 . 15 :

 315 Figure 3.15: Snapshots of simulation results of polymers of length 1000 monomers at temperature T = 0.83 and strain values of (a) 0, (b) 0.40, (c) 0.45 and (d) 0.50. The color coding is based on the local orientation order parameter S.

Figure 3 .

 3 Figure 3.16: A sequence of polymer conformations at T = 0.83 and strain-rate 10 -6 τ -1 for strains as indicated.

Figure 3 . 17 :

 317 Figure3.17: The stress-strain curves at different constant temperatures for the chain length are of 100 (a, b) and 1000 (c, d) monomers. The plots (a, c) correspond to the true strainrate of 10 -5 τ -1 and (b, d) to 10 -6 τ -1 . The points are the measured strain values and the curves are fitted using the best piece-wise linear fits, with at most three segments.

Figure 3 . 18 :

 318 Figure 3.18: The strain at which the stress changes from the initial regime of weak linear response and the start of the induced crystallization as a function of the temperature.

Figure 3 . 19 :

 319 Figure 3.19: The fraction crystallinity χ (S > 0.8) as function of simulation time at different constant temperatures for short polymer systems.

Figure 3 .

 3 Figure 3.20: (a) The radial distribution function g (r), (b) the two-dimensional perpendicular distribution g ⊥ (r), and (c) the one-dimensional parallel distribution g (r) for long chains at T = 0.86. The different strains indicated corresponds to the strain values during a constant true strain-rate of 10 -6 τ -1 . The parallel and perpendicular directions are with respect to the straining direction.

Figure 3 . 21 :

 321 Figure 3.21: The crystallinity and the average bond order P 2 along the strain direction as a function of the imposed strain (constant true strain-rate 10 -6 τ -1 ), [(a) and (c)] for short chain systems and [(b) and (d)] for long chain systems at different temperatures.

Figure 3 . 22 :

 322 Figure 3.22: The average stem length and the fraction of particles contained in stems as a function of the strain (slow strain rate of 10 -6 τ -1 ), [(a) and (c)] for short and [(b) and (d)] for long polymers at different temperatures.

Figure 3 . 23 :Figure 3 . 24 :

 323324 Figure 3.23: The number of hair-pins per polymer as function of the imposed strain at slow strain-rate, 10 -6 τ -1 at different temperatures.

Figure 3 . 25 :

 325 Figure 3.25: The (a, c) end-to-end distance R ee and (b, d) the radius of gyration parallel (R g, , solid lines) and perpendicular (R g,⊥ , dotted) to the strain direction as a function of true strain at different constant temperatures, respectively. Plots (a, b) for short polymers and (c, d) long polymers for the strain-rate 10 -6 τ -1 .

Figure 3 . 26 :

 326 Figure 3.26: The relative shape anisotropy for polymers of (a) 100 and (b) 1000 monomers during constant true strain-rate of 10 -6 τ -1 .

Figure 4 . 1 :

 41 Figure 4.1: Examples of different mapping schemes for coarse-graining of cis-polyisoprene. The circles represent the CG beads whose coordinate are defined as the center of mass (COM) position of the atoms inside the circle except in the case of (c) where the coordinate of the CG bead is located on the cross and in (d) where the mapping center of the CG bead is located at the COM of 6 atoms indicated by dotted line. The images will be referred to (a) Type-1A, (b) Type-2A, (c) Type-2B, (d) Type-2C, (e) Type-3A and (f) Type-3B, respectively.

Figure 4 . 2 :

 42 Figure 4.2: Two types of knots that are formed in two different atomic chains are depicted in (a, b) and (c, d). For easy identification of the knots, only the carbon polymer backbone is shown. All hydrogen atoms along with the methyl group carbons are not presented. The different carbons in the backbone of a monomer, CH 2 CH C CH 2 is represented by redblue-yellow-red, respectively. In (a, c) the chains are shown as tubes without any atoms, while in (b,d) the same chains are shown with the carbon atoms of the polymer backbone.

Figure 4 . 3 :

 43 Figure 4.3: Comparison of the mass densities at different temperatures for the chain length of 30 monomers along with the experimental results. The data from experiments are presented as circles and the values are taken from Ref. [173]. The measurements of a UA model used in Ref. [101] are shown as green squares and all-atom model from [150] is shown as blue square. The results from our model (AMBER potential) are the red triangles. The dotted lines are guide to the eye.

Figure 4 . 5 :

 45 Figure 4.5: The CG potentials obtained after applying STM procedure for 5 different chain lengths. The plots correspond to the models as: (a, b) Type-1A, (c, d) Type-2A and (e, f) Type-2B, (a, c) are the quartic bond potentials, (e) is a harmonic bond potential and (b) is a quartic angle potential, (d, f) are the tabulated angle potentials. The dotted line in each plot refers to the energy one k B T .

Figure 4 . 6 :

 46 Figure 4.6: The CG potentials obtained after applying STM procedure for 5 different chain lengths. The plots correspond to the models as: (a) Type-1A, (b) Type-2A and (c) Type-2B, (a, c) are the tabulated pairwise non-bonded conservative potentials. Dotted line represents the energy value which is zero.

Figure 4 . 7 :

 47 Figure 4.7: Comparison between mass density measurements using DPD simulations and reference MD values for 5 different chain lengths. The circle and triangle symbols are used to indicate reference MD and DPD measurements, respectively. The dashed lines connecting the points are provided to guide the eye. The black dotted line indicates the experimental value for long polymers after performing simple linear interpolation.

Figure 4 . 8 :

 48 Figure 4.8: The measured distributions from the DPD simulations for 5 different chain lengths are compared with a reference MD distribution. The DPD simulations are performed for the chain lengths using the CG potentials obtained from the corresponding STM process. The reference distribution corresponds to the average measurement of four independent realizations with error bar for the system with chain length of 40 monomers. The various plots correspond to different models: (a, b) Type-1A model, (c, d) Type-2A model, (e, f) Type-2B model, (a, c, e) are the bond distance distributions and (b, d, f) are the angle distributions.

Figure 4 . 9 :

 49 Figure 4.9: The measured RDF from the DPD simulations for 5 different chain lengths are compared with a reference MD distribution. The reference distribution corresponds to the average measurement of four independent realizations with error bar for the system with chain length of 40 monomers. RDF's of all monomers (excluding 1 -2, and 1 -3 bonded particles) correspond to the models: (a) Type-1A, (b) Type-2A and (c) Type-2B.

Figure 4 . 10 :

 410 Figure 4.10: Comparison between the R g and the R ee distributions obtained from the DPD simulations for 5 different chain lengths. The plots correspond to the models: (a, b) Type-1A, (c, d) Type-2A and (e, f) Type-2B, (a, c, e) are the radius of gyration and (b, d, f) are the end-to-end distributions.

Figure 4 . 11 :

 411 Figure 4.11: Comparison between the average R g from the DPD simulations and the reference MD values for 5 different chain lengths. Note that the error-bars in the MD measurement are shown only for the chain length of 40 monomers. The circle and triangle symbols are used to indicate reference MD and DPD measurements, respectively. The dashed lines connecting the points are provided to guide the eye.

3 ]Figure 4 . 12 :

 3412 Figure 4.12: Comparison between the mass densities measured using the DPD simulations for different chain lengths with nearly 100, 000 beads. The symbol (black triangles) indicates the reference MD value from atomic system with 800 monomers. The dashed lines connecting the points are provided to guide the eye. The black dotted line indicates the experimental value for long polymers after performing simple linear interpolation.

Figure 4 . 13 :

 413 Figure 4.13: Comparison between the R g distributions from the DPD simulations for different chain lengths at 350K. The distributions are scaled with the N m where N m is the number of monomer per chain .

Figure 4 . 14 :

 414 Figure 4.14: The measured MSD of COM and beads from the trajectories obtained from the DPD simulations for 5 different chain lengths are compared with their corresponding reference MD measurement. The plots correspond to the models as: (a, b) Type-1A, (c, d) to Type-2A and (e, f) to Type-2B where (a, b, c) relate to the mean square displacement of the center of mass of the polymers and (b, d, e) to the MSD of CG beads, respectively. Note that the CG simulations are performed using the CG potentials obtained for the respective chain lengths and the CG mapping schemes. The combined CG potential is used only for the chain length of 40 monomers system.

Figure 4 . 15 :

 415 Figure 4.15: Volume per monomer and average radius of gyration as function of the imposed temperature for different chain lengths. The plots correspond to the models: (a, b) Type-1A, (c, d) Type-2A, and (e, f) Type-2B, (a, c, e) are the volume per monomer and (b, d, f) are the average radius of gyration as function of the imposed temperature. Note that the CG simulations are performed using the CG potentials obtained for the respective chain lengths and the CG mapping schemes. The combined CG potential is used only for the chain length of 40 monomers system.

Figure 4 . 16 :

 416 Figure 4.16: The average volume per monomer as function of the imposed temperature for different chain lengths using the combined potential at the cooling rate of 2 K/ns.

Figure 4 . 17 :

 417 Figure 4.17: The probability distributions of (a) bond length and (b) angles for the chains length 400 beads during the constant cooling at the rate 2 K/ns in different temperature windows.

Figure 4 . 18 :

 418 Figure 4.18: The radial distribution functions of (a) all monomers, (b) the inter chain monomers, (c) the 2D perpendicular correlation and (d) the 1D parallel correlation for chain length of 400 monomers at different temperatures during the constant cooling rate 2 K/ns.

Figure 4 . 19 :

 419 Figure 4.19: The measurements of (a) the average R g and (b) fraction of bonds with (S ≥ 0.6) for different chain lengths as a function of the imposed temperature during a constant cooling at the rate 2 K/ns.

Figure 4 . 20 :

 420 Figure 4.20: Snapshots of a system with the chain length of 400 monomers during cooling at the constant rate 2 K/ns, (a) in the melt at T = 350K and (b) semi-ordered at T = 250K. The color coding is based on the local orientational order parameter S.

Figure 4 . 21 :

 421 Figure 4.21: The pressure tensor components as function of the imposed temperature during the constant cooling at the rate 2 K/ns for the chain length of 400 beads. The plots corresponds to: (a) isotropic, (b) anisotropic and (c) triclinic methods to control the pressure components.

Figure 4 .

 4 Figure 4.22: (a) The average volume per monomer and (b) the average R g as function of the imposed temperature for the chain length of 400 beads using different methods to control the pressure components at the cooling rate of 2 K/ns.

Figure 4 . 23 :

 423 Figure 4.23: Snapshots of simulation results of polymers of length 400 monomers at temperature T = 300K and strain values of (a) 0, (b) 2, (c) 4 and (d) 6. The color coding is based on the local orientation order parameter S.

Figure 4 . 24 :

 424 Figure 4.24: The stress-strain curves at constant temperatures (a) 350K and (b) 300K of different chain lengths under the true strain rate of 10 -7 /fs. The points are the measured average stress values and the curves are best piecewise continuous fits to guide the eye.

Figure 4 . 25 :

 425 Figure 4.25: The stress-strain curves at constant temperatures (a) 350K and (b) 300K of chain length of 400 monomers under different true strain rates. The points are the measured average stress values and the curves are the piecewise continuous fits to guide the eye.

  4.26.

Figure 4 .

 4 Figure 4.26: (a) Snapshot of a configuration with the primitive path and (b) a single polymer with its corresponding primitive path in red color where the neighboring polymers are hidden. Note, the presented images correspond to the chain length of 400 monomers.

Figure 4 .

 4 Figure 4.27: (a) Kuhn length and (b) entanglement length for different chain lengths at two temperatures. The dashed lines connecting the points are provided to guide the eye. The experimental values for long polymers are obtained after performing simple linear interpolation.

Figure 4 . 28 :

 428 Figure 4.28: The entanglement lengths as function of the imposed strain at constant temperature (a) 350K and (b) 300K of different chain lengths under the true strain rate of 10 -7 /fs.

Figure 4 . 29 :

 429 Figure 4.29: The probability distributions of (a) bond length and (b) angles for the chains length 400 beads under the true strain rate of 10 -7 /fs measured in different strain windows.

Figure 4 .

 4 Figure 4.30: (a) The inter-chain radial distribution function g(r), (b) the 2D perpendicular distribution g ⊥ (r), and (c) the 1D parallel distribution g (r) for the long chains of 400 monomers at T= 300K and different strains during constant straining rate 10 -7 /fs. The perpendicular and parallel directions with respect to the strain direction.

Figure 4 . 31 :

 431 Figure 4.31: The fraction number of bonds with (S ≥ 0.6) as function of strain at constant temperature (a) 350K and (b) 300K of different chain lengths under the true strain rate of 10 -7 /fs.

Figure 4 . 32 :

 432 Figure 4.32: The local orientational order distribution as function of strain at constant temperature 300K of chain length 400 beads per polymer under the true strain rate of (a) 10 -6 /fs, (b) 10 -7 /fs and (c) 10 -8 /fs. Each distribution represents a strain window with 20 configurations during straining.

Figure 4 . 33 :

 433 Figure 4.33: The average bond order P 2 along the strain direction as a function of strain at constant temperature (a) 350K and (b) 300K of different chain lengths under the true strain rate of 10 -7 /fs.

Figure 4 . 34 :

 434 Figure 4.34: The radius of gyration along (R ,g , solid lines) and perpendicular (R ⊥,g , dotted lines) to the strain direction as function of strain at constant temperature (a) 350K and (b) 300K of different chain lengths under the true strain rate of 10 -7 /fs.

Figure 5 . 2 :

 52 Figure 5.2: Mass densities for different chain lengths with different cross-link densities measured during DPD simulations. The black triangle symbols are used to indicate reference MD measurements. The dashed lines connecting the points are provided to guide the eye. The black dotted line indicates estimated experimental value for long polymers after performing simple linear interpolation based on existing data.

Figure 5 . 3 :

 53 Figure 5.3: The bond length distributions of chain length 400 beads under the true strain rate of 10 -7 /fs with (a) no cross-links, (b) 1% cross-links and (c) 2% cross-links measured in different strain windows. The vertical dotted line corresponds to the bond distance of 5.5 Å.

Figure 5 . 4 :

 54 Figure 5.4: The stress-strain curves at constant temperature 350K for (a) different chain lengths with 1% cross-link density and (b) chain length of 400 monomer system with different cross-link densities under the true strain rate of 10 -7 /fs. The points are the measured average stress values and the solid line presented for all the curves is guide for the eye.

Fig. 5

 5 Fig.5.4b depicts the stress-strain plots for the 400 monomer chain system with varying cross-link densities. It can be noticed that for cross-link density ≥ 2%, the strain range is truncated based on the bond length cutoff criteria. It can be observed that the higher the cross-link density smaller is the strain range where the linear stress response occurs. This initial domain is followed by the divergence of the stress, i.e., the non-linear stress-strain behavior. The divergence of stress after reaching a particular strain is observed in many experimental works for vulcanized IR[177,182, 190]. When the system has more crosslinks, a path between a monomer and its image over the periodic boundary condition will exist, i.e., a percolating network is formed. For higher cross-link density systems, the system realizes it is a percolating network at smaller strain this results in the shorter initial regime where the stress response is linear.To analyze the local structure at monomer scale we compute the local orientational order. We use the nematic order tensor Q αβ expressed in Eq. 3.4 and calculate the largest

Figure 5 . 5 :

 55 Figure 5.5: The fraction of bonds with S ≥ 0.6 as function of the imposed strain at a constant temperature 350K for (a) different chain lengths with 1% cross-link density and (b) chain length of 400 monomer system with different cross-link densities under the true strain rate of 10 -7 /fs.

Figure 5 . 6 :

 56 Figure5.6: The average bond order P 2 along the strain direction as function of the imposed strain at a constant temperature 350K for (a) different chain lengths with 1% cross-link density and (b) chain length of 400 monomer system with different cross-link densities under the true strain rate of 10 -7 /fs.

Figure 5 . 7 :

 57 Figure 5.7: The radius of gyration along (R ,g , solid lines) and perpendicular (R ⊥,g , dotted lines) to the strain direction as function of the imposed strain at a constant temperature 350K for (a) different chain lengths with 1% cross-link density and (b) chain length of 400 monomer system with different cross-link densities under the true strain rate of 10 -7 /fs.

Figure 5 . 8 :

 58 Figure 5.8: The end-to-end distance R ee as function of the imposed strain at a constant temperature 350K for (a) different chain lengths with 1% cross-link density and (b) chain length of 400 monomer system with different cross-link densities under the true strain rate of 10 -7 /fs.

  .19. It can be seen that the crystallinity increases for the systems at T < 0.82 over time resulting in a semi-crystalline phase without applying any mechanical loading.

	0.78	0.8	0.82	0.84	0.86	0.88	0.9	0.92

Table 4 .

 4 1: The measured densities for atomic systems for five different chains lengths at temperature 350K and atmospheric pressure.

	Chain length [monomers]	ρ [g/cm 3 ]
	10	0.840
	20	0.855
	30	0.861
	40	0.863 ± 0.001
	80	0.867

Table 4 .

 4 2: Computed mean square radius of gyration and mean square end-to-end distance for five different chain lengths at temperature 350K and atmospheric pressure.

	Chain length [monomers] 〈R 2 g 〉 [Å 2 ]	〈R 2 ee 〉 [Å 2 ]	〈R 2 ee 〉 / 〈R 2 g 〉
	10	49	303	6.2
	20	108	666	6.2
	30	150	932	6.2
	40	187 ± 20 1040 ± 200	5.6
	80	415	2447	5.9

Table 4 .

 4 3: The number of interaction parameters for different CG models provided to STM for optimization.

	CG model	Bond	Angle	Pairwise	Total number
	types	parameters parameters parameters of parameters
	Type-1A	4	4	20	26
	Type-2A	4	33	20	55
	Type-2B	2	33	20	53

Table

Unpublished

was carried in a joint public-private laboratory, Simatlab, between the Clermont Auvergne University and the company Michelin. First of all, I would like to thank the Clermont Auvergne Project (CAP) 20-25 for funding this thesis and all the members of Simatlab and l'Institut de Chimie de Clermont-Ferrand who have helped me complete this thesis.

carbon-carbon double bonds. As the carbon-carbon single bond is included within a CG bead, the different conformations due to the rotation around the single bond is suppressed resulting in a skewed uni-modal distribution. We use a class2 bond potential to model the bond interaction similar to Type-1A option. The angle distribution for this mapping choice has a double peak profile with the first peak close to 90 • and another peak near 140 • .

To obtain a CG model for this mapping scheme, the angle interaction is modeled using a tabulated style.

The measured structural properties using the reference trajectory after applying the Type-2B mapping scheme have uni-modal and bi-modal profiles for bond length and angle distributions, respectively. Since the mapping center of the CG bead is located on the carbon-carbon single bond, it does not depend on the dihedral angle of the 4 carbon atoms along the monomer backbone. This results in the bond distance profile having a clear single the DPD simulations have wider distributions for all the CG models. It can be noticed that all the chain lengths in the DPD measurements for Type-1A and 2A are visually alike. Here, the location of the CG beads are at the center-of-mass for the groups of particles considered within their respective CG beads. This results in a wider bond length distribution as shown in Fig. 4.8a,c for Type-1A and 2A, respectively. Due to this, the effect of the end monomers on the short chains are not visible.

In the Type-2B model the CG bead is located at the center of the carbon-carbon single bond. Due to this, the bond length is clearly defined as compared to the other CG mapping schemes. This results in the bond length distribution having clear single-peak as shown in Fig. 4.8e. For the short chains of 10 monomers the number of middle monomers is relatively smaller than for the longer chains considered here. The combination of clearly defined bond distance and the effect of end monomers, manifests in wider bond length distribution for the chain length of 10 monomers. One possible way to confirm the hypothesis of the effect of end monomers is to measure the bond length distributions of end-and middlebonds separately. If we were to measure such distributions, we expect that the bond length distribution for end-and middle-bonds will be same for both short and long chains but the bond length distribution that include the contributions of both end-and middle-bonds will be chain length dependent.

Both the quartic (Type-1A) and the tabulated (Type-2A) angle potentials do well in reproducing the reference angle distributions. In the Type-2A model, the measured angle distribution for the long chains appear closer to the reference profile, but the distributions of other chain lengths are not very different from the 80 monomers system. In the Type-2A and 2B case, the angle distributions have a bi-modal shape. The nature of the curves irrespective of the chain lengths appear to be similar to the reference distributions with more population at the first peak than the second one. In the Type-2B model, different chain lengths have different probabilities in the DPD case, but a clear chain length dependency cannot be established just using the data presented here.

Overall STM does a reasonably good job in reproducing both angle and bond length distributions given the fact that the Bayesian method does not strictly focus on the structural measurements unlike other CG methods such as IBM or IMC.

The short-range structures is investigated by computing the radial distribution function for all monomers excluding 1 -2, and 1 -3 bonded particles. The RDF's for different chain lengths obtained from the DPD simulations and the average MD measurement of chain length 40 monomers obtained from four independent realizations for different CG models are shown in Fig. 4.9.

The location of peaks are well reproduced in the DPD measurements when compared to the reference RDF in Type-2B model. In the Type-1A and 2A models the RDFs do not reproduce the second peak in the DPD measurement for the chain length of 40 monomers compared to the reference MD shown in Fig. 4.9. This observation holds true for other chain lengths considered for which the reference MD measurements are not presented here.

In the Type-1A model, the density is in good agreement along with the bond length and the angle distributions compared to the MD measurements. The reference bond distance distribution is bi-modal, which means if we use this model and perform cooling or straining simulation, the bond distances can easily change from short to long distance. The consequence of preferred longer bonds is that upon cooling the volume of the system will increase which is undesirable. Hence, Type-1A model cannot be used to study crystallization of cis-PI The radius of gyration and the end-to-end distance measured at two temperatures are similar for different chain lengths. For the short polymers of length 40 monomers the average end-to-end distance and the average R g component (parallel to the straining direction) increases slightly. Here the system remains a melt and the mobility is large enough to eliminate induced deformation in the overall polymer shape. In the case of long chains of 400 monomers, the average R ,g increases from nearly 30 Å to about 180 Å at maximum strain for both temperatures. Similarly, the end-to-end distance at zero strain is about 110 Å and at the maximum strain the average value increases by a factor 3.5. On closer observation we can notice that the slope of this curve changes approximately at the same strain value at which the local orientational order starts to increase. For the long polymers, the timescale required to untangle exceeds the simulation time and hence the polymers are forced to respond to the imposed strain. To illustrate the effect of the strain rate, i.e., the simulation time, we present the 〈R g 〉 and the 〈R ee 〉 as function of the imposed strain for the chain length of 400 monomers in Fig. 4.36.

Here we can see both the measurements increase as function of strain with the fast strain rate resulting in higher average value and slowest strain strain rate in the smallest average value.

Summary and discussion

We have performed atomic simulations of cis-PI model of different chain lengths at 350K. The results from the atomic simulations were compared with other atomic models [101,150] and experiments [173] from literature. We find that our atomic model is capable of reproducing the density within 1 % deviation from the measurement found in literature [101] at temperature 350K and 413K. Based on these results we can say that our atomic model is good to reproduce the densities with respect to the mentioned literature studies.

In the next step, we apply different CG mapping schemes to the atomic trajectories. These geometrically mapped reference trajectories were subjected to the STM process to obtain various CG potentials for different chain lengths. Using these CG potentials, DPD simulations were performed and various structural measurements were computed. To validate out CG models, we compare our DPD results with atomic measurements. The mass densities measured using CG simulations are within 1 %, 4 %, and 3 % deviations for Type-1A, Type-2A and Type-2B CG models for chain lengths ≥ 30 monomers, respectively. Further, the measured bond length, angle and RDF (excluding 1-2 and 1-3 bonded particles) from the CG simulations are compared with respect to the distributions measured from geometrically mapped reference trajectories. From the results concerning the bond length and the angle distributions we can say that the STM potentials are capable of reproducing the structural measurements (bond lengths and angles) well when compared with the reference distributions given the fact that the Bayesian method does not focus on the structural measurements.

All the CG models that we have considered in our study are simple approaches to model cis-PI and therefore, each CG model has its own limitations.

We recall that, although the density and other structural measurements like bond length, and angle distributions are well reproduced in the case of Type-1A model the bond distance distribution is bi-modal. If we use this CG model to perform cooling or straining simulations, the bond distances can easily change from one location to other. Hence, this model cannot be used to study crystallization. In Type-2A CG model bond has a skewed profile with one peak and angle distribution is bi-modal. These distributions are well reproduced, however, all the peaks of g(r) are not well reproduced. When we tested this model we did not see any sign of phase transition. Based on the g(r) (excluding 1 -2 and 1 -3 bonded particles) we can say that the CG potentials obtained from STM reproduces the location of peaks well in the Type-2B CG model compared to the Type-2A CG model.

We again state that in the Type-2 variants (A, B) the end atoms are not considered. One

Appendix A

Interaction potentials A.1 Atomic force-field parameters: AMBER

The force-field parameters used to perform atomic simulations of cis-PI system are provided in this section.

A.2 CG potentials

The CG potentials obtained from STM fitting procedure is presented below. This is a combined potential obtained from four independent realizations of the system with 40 monomers per chain. We have used harmonic bond potential with the parameters: Cis-Cis N 20