An atomistic molecular dynamic model to study the properties of LLDPE and wax

Wax is often physically mixed with linear low-density polyethylene (LLDPE) to form a new polymer material. However, the morphology of these materials has not been described sufficiently, because molecular interaction at an atomic level was inadequately studied. Molecular dynamics (MD) simulation, using the Material Studio software as a computational tool, was available to develop models for wax and LLDPE to study their properties at an atomistic level. The models were validated by comparing the properties, such as solubility, density, and transition temperature, obtained with the models, with those obtained from experiments. After validation, the application of the model showed that the branch content of modelled LLDPE affected the glass transition temperature when the branch content was 70 per 1000 carbons for models with different branch lengths. The longer the branch length, the higher the glass transition temperature of LLDPE. However, the solubility parameter was unsuccessful in finding the length of LLDPE required to represent a single chain, because the properties of a copolymer are affected by the length of the chain, the number of branches, and their distribution on the backbone. The chain length of the wax showed no relationship to the solubility parameter in the solid state or in the melt. There was a decrease in the solubility parameter of the modelled LLDPE with an increase in temperature. The LLDPE and wax properties, examined through MD simulations, were within 10% of the experimental values.


Introduction
In a continuous search for better polymeric materials in terms of cost and tailor-made properties, wide-ranging applications and easy to recycle, it was found that polymer blends and composites are the solution [1,2].The physical mixing of different polymers and wax or several polymers is regarded as the cheapest method of creating novel polymeric materials with specific properties not possessed by the individual polymers or the wax [3,4].
Although Luyt and co-workers [5][6][7] studied many properties of various polyethylene/wax blends, limited success has been achieved in establishing the exact morphologies responsible for these properties.For example, experimental studies showed that wax is partially miscible with LLDPE at low content; however, the interaction at a molecular level is unknown [8,9].Therefore, the exact morphology of the LLDPE/wax blends has not yet been fully described.This lack of morphology information is probably because the wax' structural and physical appearance is almost exactly like that of LLDPE.Due to the similarity in appearance, it is difficult for both transmission electron microscopy (TEM) and scanning electron microscopy (SEM) to show the exact arrangement or interaction between the wax and LLDPE components on the surface or inside the blend.
Understanding the miscibility of the polymers or polymer and wax in the polymer blend is the first step in establishing the exact morphologies and predicting the physical properties of polymer blends [4,10].Polymer blends are grouped into three groups, namely immiscible, compatible and miscible.Miscible blends have one continuous phase and one 416 Page 2 of 16 glass transition temperature, while immiscible blends have more than one glass transition temperature, depending on the length scale and number of components in the blends [2,10].Therefore, one way to determine the miscibility of blends is to investigate the glass transitions of the pure polymers, and to compare them to the glass transition temperature(s) of the blends.
The compatibility in polymer blends could be investigated by the Flory-Huggins interaction parameter (χ) [11,12], which measures the interaction between polymers defined by the excess in free energy of mixing of the respective polymers.The Flory-Huggins interaction parameter was determined through various experiments, such as inverse gas chromatography [13], heats of mixing [14,15], comparisons of solubility parameters [16], and computational simulations using force fields in molecular dynamics [17].The thermodynamic properties of blends, such as the energies of mixing, polymers' densities, glass transition temperatures, and solubility parameters, are calculated by using molecular dynamics simulations [18][19][20].Many molecular simulation studies were done to predetermine the blends' morphology properties by simulating the compatibility between the components of several binary blends [21,22].However, before molecular dynamics can be used to study the blends, models of the pure components and the blends must be developed and validated, as molecular structures are formed by atoms.
In this paper we report on the use of molecular dynamics simulations to develop and validate models of wax and LLDPE through comparison of the models' properties, such as solubility, density, and transition temperature with that of the same materials' experimental properties.

Experimental
This paper covers the force field based molecular simulation of pure LLDPE [23] and synthetic paraffin wax, obtained via the Fischer-Tropsch process [24], for model development and validation.The model construction and development involved the creation of the molecular structures of the two monomers that form LLDPE, namely ethylene and 1-hexene.The created monomers and the LLDPE molecules were further optimised by the DMol 3 module [25] in the Materials Studio software [26] to relax the structures to the lowest minimum energy.Different lengths of LLDPE with different branch contents were constructed using these monomers.Using the Amorphous Cell module [26], one or more LLDPE molecules were added into three-dimensional periodic systems (unit cells).The resulting systems were then further relaxed using the canonical constant number, volume and temperature (NVT) and the isothermal-isobaric (constant number, pressure and temperature) (NPT) ensembles at 298 K and 500 K, followed by annealing the unit cells from 100 to 500 K and back to 100 K to further equilibrate the systems.
Properties such as solubility (through calculation of the solubility parameter), density at various temperatures, radius of gyration, and radial distribution function were evaluated from the final 100 dynamics (NPT) steps.Density vs temperature graphs were used to determine the glass transition temperature of the modelled LLDPE with different branch lengths and branch content [27].

Molecular simulation procedure
Materials Studio (MS) 2019 modelling software (BIOVIA © ) [26] made available by the Centre for High Performance Computing (CHPC) [28] situated at the Rosebank campus of the CSIR (Council for Scientific and Industrial Research) in Cape Town, South Africa, was used for the modelling calculations.The CHPC supercomputer cluster called Lengau (1368 compute nodes, 5 large memory nodes, 9 GPU nodes, 4 PF storage over the Lustre parallel file system) was used to run all the molecular simulations.Between 48 and 96 cores of the CPU nodes were required to run in parallel on the clusters for a sufficient simulation time specifically to do molecular dynamics.A personal computer (PC) (Dell Opti-Plex 9020, 2xQuad core Intel i7-4770 s 3.10 GHz, 16 GB RAM, 1TR GB HDD), situated at the University of the Free State (Qwaqwa Campus), was utilized as the workstation to upload calculations, and receive the downloading of the results from the cluster.

Construction and optimisation of monomer models
The monomers of LLDPE, ethylene and 1-hexene (Fig. 1) were drawn manually using the Visualiser module of the Materials Studio (MS) software [26].The two monomers listed here are used as the starting compounds for the synthetic catalytic polymerisation of LLDPE under specific pressure and temperature conditions.Since the polymerisation of LLDPE occurs according to a free radical mechanism [29] with the addition of an organic peroxide radical initiator, the two monomers were optimised as radical molecules by changing the multiplicity to a triplet to simulate one unpaired electron to allow for the formation of the polymers in the MS software.
The two monomers were further relaxed by density functional theory (DFT) optimisation as implemented in the DMol 3 module [30][31][32] of the Materials Studio software [26].Geometric optimisation was conducted using the generalised gradient approximation functional of Perdew and Wang (GGA-PW91) [33] incorporating a double numeric Page 3 of 16 416 polarised (DNP+) and an additional diffuse function basis set.This basis set was satisfactorily tested against polymers and surfaces [34], especially for calculations requiring a large orbital cut-off [30].Additionally, less simulation time is required when DNP+ basis set is opted for especially simulations which pre-require experimental validations [35].
The optimisation process of LLDPE and its two monomers' structures was carried in the following set of steps: medium was chosen as the optimisation quality, at the energy threshold of 2.0 × 10 -5 Hartree (in the Hartree wave function) maximum and force of 0.004 Hartree Å −1 , with a possible 50 iterations in total, the length between Cartesian coordinates (maximum allowed length in search of Cartesian coordinates) was at a size of 0.3 Å, and between points 0.005 Å displacement.Because of the unpaired electron belonging to the radical monomers, triplet multiplicity state was apparent, and symmetry was unselected.The following were the electronic setup settings: integration accuracy was put to medium, 1.0 × 10 -5 tolerance was opted for selfconsistent field (SCF), with 50 SCF cycles maximum, as well as multipolar expansion (hexadecapole) with a global orbital cutoff set to medium and smearing of 0.005, and all electron core treatment was chosen.The optimisation calculations of the structural molecules used were done in the absence of solvents.The only properties determined from the monomer structures after optimisation were the energy evolution and convergence.

Construction and optimisation of LLDPE using a force field
The two optimised monomers were used to create different single LLDPE chains with varying lengths and branch content (Fig. 2).The random copolymer interface, which is found in the built polymers of the MS software [26], was used to create the LLDPE chains.However, the polymerise option within a random copolymer interface made it possible to first insert the two constructed monomers into two separate columns of the template.The polymerise tab has many other options, such as using forced concentrations or reactivity ratios (which allows variation of branches along the polymer backbone) and chain length and/or the number of chains that determine the size and shape of the polymer molecule.This construction was done for chain lengths of 10 to 500 [36] with 1-hexene as branches (varying content) distributed randomly along the polymer chain backbone [36,37].These LLDPE polymer chains were geometrically optimised within the Forcite module by the force field called condensed phase optimised molecular potentials for atomistic simulation studies (COMPASS) [38].Geometry optimisation was run through a Quasi-Newton algorithm in order to quickly reach the global minimum (Fig. 3).The convergence tolerance setup was as follows: high quality (fine) optimisation for short LLDPE chain structures, with the energy at 2.0 × 10 -5 kcal mol −1 and maximum force of 2.0 × 10 -5 kcal mol −1 A −1 , different maximum 50 iterations and 0.005 Å displacement.The long chains (Fig. 2b) were optimised at an energy tolerance of 1.0 × 10 -3 kcal mol −1 and 0.5 kcal mol −1 A −1 of force, and 90,000 maximum number of iterations.The Gasteiger charge equilibration method was used to set up parameters for the charge contribution by all the atoms in the LLDPE model during the optimisation process.The maximum number of iterations per cycle was 50 for the charge calculation.The convergence limit, which specifies the total transfer of charge permitted between two atoms apart, was 2.0 × 10 -5 .The original publication on the Gasteiger method [39] covers parameters for the halogens and the common elements such as carbon, hydrogen, oxygen, nitrogen, and sulphur.

Construction and optimisation of wax
The MS software contains the readily available ethylene molecule which was used to produce the wax chains and not LLDPE, which required ethylene to be drawn manually.The created wax chains of varying lengths (C 20 H 22 -C 44 H 46 ) were optimised by COMPASS [38].The chain length of the wax was limited to C 20 H 22 -C 44 H 46 because soft medium wax was used in the experiments [36].The wax chains geometry optimisation was carried through quasi-Newton algorithm and the energy convergence was quite quick, especially for short polyethylenes.The Amorphous Cell module incorporated in MS was quite useful for the creation of the amorphous cells in which to put wax chains for further MD simulations.These cells containing the wax chains (periodic systems) were reduced in size by the canonical constant number, volume and temperature (NVT) ensemble MD simulations, first at 298 K for a total simulation of 1000 ps and 1 fs time step, then at 500 K.The resulting cells were submitted to an isothermal-isobaric simulation (constant number, pressure and temperature) (NPT), at 298 K and 1 fs time step for a total simulation of 1000 ps, and then a final NPT simulation using the same settings as in previous step was conducted at 500 K.

Building three-dimensional polymeric systems using the amorphous cell module
Amorphous cells of LLDPE molecules were also produced in order to further minimize the chains before MD simulations were started in Materials Studio.The LLDPE and wax models were put inside the amorphous cells by the common technique detailed by Theodorou and Suter [40].The calculation option in the Amorphous Cell module takes a short time during construction, and it was preferred over the construction (legacy) option that has been used in the older versions of the software.The calculation option settings followed these steps: cubic was opted as lattice type, ramping the initial density of 0.600 g cm −3 to the final densities of 0.900 g cm −3 for wax and 0.924 g cm −3 for LLDPE.These density values were chosen because they are densities of the neat polymers prepared experimentally [36].During the construction, torsion angles were determined while relaxing the chains without cell optimisation.An energy convergence level was set to 2.0 × 10 -5 kcal mol −1 for smaller molecules and 1.0 × 10 -3 kcal mol −1 for large molecules through quasi-Newton.The van der Waals was preferred for all energy calculations while the Ewald was used for electrostatic summation, the original method for calculating charges was Gasteiger (Gast_original1.0),and COMPASS was selected over other types of force fields.

Molecular dynamics (MD) equilibration
One molecule chain per cubic simulation box (periodic system) was constructed in the case of long LLDPE chains (10 to 100 branch content for a fixed 1000 carbon atoms in the LLDPE chain) as shown in Fig. 2. The molecule inside the periodic systems was relaxed through initial canonical constant number, volume and temperature (NVT) and constant number, pressure and temperature (NPT) MD ensembles to equilibrate the periodic system.Both NVT and NPT relation dynamics were done first at 298 K and 1.01 × 10 -5 GPa with a 1.0 fs time step, using 4.0 × 10 5 steps for a total simulation of 1000 ps.A random initial velocity was applied in the whole simulation, using the non-linear Hebbian learning (NHL) algorithm and Berendsen barostat to keep the volume or pressure constant.
A similar process was followed for the case where a cubic simulation box contained more than one molecule, to look at the effect of the molecules content per amorphous cell relative to the total density.This process was followed for short chains of LLDPE (whereby the branch content was 1 to 25 for 100 carbons long) and wax.The MD simulation calculations were created in the absence of a solvent.
The chains of both wax and LLDPE in the cells were minimized by an annealing simulation procedure from 150 to 450 K and back to 150 K with a heating ramp of 0.1 K ps −1 for 30 cycles.The minimisation annealing started with NVT, then followed by NPT MD [20,41,42] while the time step was still set at 1.0 for a total simulation of 1.0 × 10 5 steps.Torsion angles were predetermined at 0.25 radii of Van der Waals forces for a 1000 ps loading steps.The quasi-Newton method was chosen as the algorithm method for 1000 ps with 50 iterations maximum and convergence limit of 5.0 × 10 -6 .The charge equilibration method was implemented to evaluate the partial charges of the atoms involved.Electrostatic energy between the atoms in the cells was considered by the Ewald summation method at a limit of 1.0 × 10 -4 kcal mol −1 (fine).The van der Waals interaction was selected for the atom-based summation with a cut-off distance of 15.5 Å.When the 30-cycle annealing process was completed, further 1000 ps NVT MD simulations were carried out for 4.0 × 10 5 steps as the final relaxations of the amorphous cells at 1.0 fs step time, followed by NPT MD at 298 K first then finally at 500 K.
From cubic amorphous cells, supercells (unit cells extending three times in all three dimensions) were constructed with a different number of polymer or wax molecules.The created supercells were relaxed further with the NVT MD simulations for short LLDPE chains and wax models to create new cells with 27 sub-units.These supercells were first minimised with NVT at 298 K, and pressure of 1.01 × 10 -5 GPa and a total simulation of 1000 ps at 1.0 fs step time.NPT dynamics were used under the same conditions to find the models' densities, and then the models were further relaxed.A random initial velocity was applied in all the simulations, and a 0.01 q ratio was used through an NHL algorithm and Berendsen barostat.
A production procedure was constructed by molecular dynamics on the LLDPE single cells, or supercells using an NPT ensemble from 500 to 100 K at the total simulation of 1000 ps and a 1.0 fs step time without changing any settings to simulate the slow cooling of the pure modelled LLDPE.Glass transition temperature (T g ) was determined from the graphs generated when various thermodynamic properties were obtained from the results.The three properties which were obtainable to determine T g were the changes in density, volume, and thermal expansion versus temperature [27,42].
Density was again used with the Hildebrand solubility parameter (δ) through cohesive energy density calculation to validate the LLDPE models.Solubility of the LLDPE as a function of branch content and branch length was determined by calculating cohesive density after molecular dynamics.Density versus temperature was used to investigate the thermal behaviour of LLDPE and wax [43].

Geometry optimisation of structures
Geometric optimisation is the starting point for any atomistic molecular dynamics study, whether it is a small or a large system.Geometric optimisation aims to get the structures to the lowest possible energy configuration [44].This optimisation is a continuous movement of atoms until the forces between them are zero or at a minimum.
Many algorithms have been developed to minimise small molecules through quantum mechanical calculations, while others are suitable for molecular mechanics' minimisation of larger molecules [45].In molecular mechanics, optimisation is done on large systems, and the potential energy function will have too many variables to solve.Therefore, energy minimisation was undertaken using density functional theory (DFT), since it predicts more accurately the electron density of the molecules, especially if the information on how materials are structured and behave is required [46].Energy minimisation can also be solved by molecular mechanics using a force field if large systems are modelled or a high performing computer with a large central processing unit (CPU) is required.Forcite is such a force field method that could be used to minimise structures, while DFT calculations are done using the DMol 3 module to optimise structures when using Materials Studio software [26].
The two monomers (ethylene and 1-hexene) used for building LLDPE polymer structures were optimised using DMol 3 .It was expected that 1-hexene would be minimised slower than ethylene because more atoms lead to a significant increasing number of parameters, all of which must be minimised (Fig. 3).In molecular mechanics, energy minimisation of molecules is usually achieved by using different types of algorithms.The nature of the problem at hand dictates the type of algorithm that must be used.The energy minimisation algorithm methods are divided into the first and second derivatives [47].
There are five different algorithms available in the Materials Studio software for the Forcite module.It was therefore imperative to explore the most suitable algorithm that works best for both wax and LLDPE structures, because of the sizes of the systems that would have to be used to carry out simulations on pure structures and blends.It was therefore important to consider factors such as step size or time to carry out the simulations in order to do the simulation with a small chance of calculation failure.Using all five algorithms, the best algorithm to optimise the LLDPE model with 20 carbons in the backbone was identified by optimising LLDPE under the same settings (Fig. 4).
A closer look at Fig. 4 shows that the energy versus optimisation steps go to a minimum more rapidly when the conjugate, adopted basis Newton-Raphson (ABNR), and quasi-Newton methods were applied.The energy was minimised to -30 kcal mol −1 within the first 300 optimisation steps for all three algorithms.The conjugate method corrects the failure of the steepest descent, and it accumulates the information about the function from the previous iteration to the next [44].It is quicker than the steepest descent, even though it Page 7 of 16 416 is a first derivative technique.The quasi-Newton and ABNR methods use the second derivative together with the gradient to search for the global minimum.However, the ABNR efficiency increased when convergence was approached and required huge computational storage or cost (Fig. 4).Looking at Fig. 5, quasi-Newton looks like it converged and minimised LLDPE within 600 steps.Because quasi-Newton seeks to build an approximation to the Hessian inverse, it converges well and does not require a line minimiser like steepest descent, and therefore it is more efficient to use [44].
The log of the convergence in both energy and gradient (as depicted in Fig. 5) of all the different algorithms shows that both Smart and Quasi-Newton minimise the energy of the LLPDE models very quickly and in short steps (small gradient).In contrast, the steepest descent algorithm was a slow but robust method that took long to converge close to the minimum.The gradient fluctuated close to the same energy value after 2000 optimisation steps (Fig. 5(b)) until it converged in 5000 steps.ABNR and the conjugate gradient took 3500 and 4000 steps, respectively, to get to the minimum energy value and with a large gradient, which means that more CPU resources were required to minimise the configurations of LLDPE.
The current state and size of the structure dictates the choice of algorithm to use for optimisation.The conjugate gradient and steepest descent are widely used methods, provided the molecules are not too far from equilibrium.The conjugate gradient was chosen because only a few previous gradients were found compared to the Newton-Raphson methods which required directions to be stored.Most Newton-Raphson methods require large disk space to store a second derivative of the matrix, therefore it was not ideal to apply them.Because each path reverses the progress made by the earlier iterations, the steepest descents converge slowly close to a minimum.Furthermore, conjugate gradient always prevents the next vector direction from undoing the earlier progress.It produced a complete set of conjugate directions in such a manner that the successive step keeps on redefining a new direction towards the minimum.Therefore, although the steepest descent was robust, it is not ideal for optimising the large polymer structures of LLDPE due to computational costs.It is the method that generates a stable minimumenergy molecule, regardless of the beginning of a process or function.It is often used when dealing with large molecules.It simply means it requires more time to get to the minimum energy in search of the local minimum when the last three algorithms are used.
The quasi-Newton-Raphson follows the basic method of conjugate-gradients.The method could successfully optimise all the LLDPE models in the shortest simulation time, because the optimisation gradients follow a Newtonian framework.

Evaluation of equilibrium structures
The minimum energy structures are usually confirmed using MD plotting of either the different energies or the density at 298 K versus the simulation time [48,49].Pressure and different energies of a wax-18 model (wax model with 18 carbons atoms) versus simulation time, as shown in Fig. 6, confirm that the wax-18 structures reached equilibrium in 250 ps.However, when the same model was systemically cooled from 350 K (above the wax melting temperature), it took at least one nanosecond for the structures to equilibrate.Therefore, the wax models equilibrated well within 1000 ps.

Validation of LLDPE and wax models
There should be some form of comparison with experimental or previously reported literature data to verify any computational study's accuracy level.This comparison is usually done by comparing the models' calculated properties against experimental measurements or results of other researchers' computational studies.The models developed for neat LLDPE and wax were validated by comparing the thermodynamic properties such as density, volume, temperature and enthalpy against simulation time, branch content and chain length.
The density of modelled long LLDPE was plotted against the number of branches per 1000 carbons in the backbone (Fig. 7(a)).There was an observable increase in density as the number of branches per 1000 carbons increased up to 50 branches, which is expected from a less branched polyethylene.Their chains tend to fold in an ordered way as a densely packed network is created.An increase in the number of branches further obstructs chains from easy folding since branching keeps individual chains further apart, reducing the material's density.Therefore, a sudden decrease in density when 50 branches or more were added to the backbone is expected.
A very linear PE chain will easily approach other PE chains that are similar in structural arrangement and create a packed network.This type of packing is the reason why a high-density polyethylene (HDPE) is relatively strong and stiff, and why LLDPE is more flexible than HDPE even though it is stiffer and stronger than low-density polyethylene (LDPE) (with long branches).The network does not say anything regarding the length of the individual molecules, or even of the branched chains and they may be short or quite long [50].If the chains are long, the melt index, will be low and the material will reveal higher melt viscosity and strength.A moulded LLDPE molecular arrangement will have enhanced impact performance than a material of the same density made up of shorter chains.The density of short chains (100 carbons in the backbone) of modelled LLDPE with different branch content is shown in Fig. 7(b).The density of LLDPE increases as the branch content of the models increases.A similar trend is seen when wax density is plotted against chain length (Fig. 7(c)).In Fig. 7(b), where the chains were shorter (100 carbons) than the longer chains (1000 carbons) of Fig. 7(a), the density increased continuously beyond the 4:100 ratio, which is the same as 40:1000 in the case of the long chains.This increase clearly shows that short chains play a different role compared to long chains in polyethylenes.This argument can be supplemented by the same observations when the density of soft paraffin wax is plotted versus chain length.
The density increased faster with the short than with the long molecular chains (Fig. 7(c)).Determination of glass transition temperature (T g ) is another way to verify that the models' properties are the same as those of the neat polymers.Therefore, the T g of both LLDPE and wax was determined using three possible methods proposed by other researchers [51,52].Changes in density and volume versus a change in temperature of the slow cooling process modelled with NPT dynamics were used to determine T g .
The variation of density versus temperature was the first approach used to investigate the T g of LLDPE with different branch content per fixed length of the chain in the backbone, as illustrated in Figs. 8 and 9.Ten to 100 branches attached to 1000 carbons long making the backbone were simulated as a function of temperature to determine the effect of branches on the simulated transition (Fig. 8).The density of a short-chain LLDPE (with 1 to 25 branches per 100 carbons) was also simulated as a function of temperature when cooled from the melt (Fig. 9).The number of branches did not show any significant difference between different chains in terms of length or branch content in the obtained models.
Figure 10 shows the thermal expansion coefficient against temperature in a systematic slow cooling process of NPT dynamic simulation.The amorphous cells of the models were cooled slowly from 500 to 100 K at a cooling rate of 25 K ns −1 per simulation time of 1000 ps.The thermal expansion coefficient of the modelled LLDPE decreased until it reached a transition at about 175 K.The same is observed in Fig. 11 when cell volume versus temperature was used.Therefore, the transition of LLDPE was the same whether density, volume or thermal expansion was used from the simulated models.
Table 1 shows the T g of several LLDPE models in the 160 to 180 K region close to the region observed in experiments [53].DMA experiments were done on LLDPE with butyl and hexyl branches (1-octene) by Simanke et al. [53,54] (Table 1).Mathews et al. [55] did DMA on butene-based LLDPE samples cooled by quenching and constant room temperate.The values of the glass transition temperature of LLDPE with less than 50 hexyl branches per 1000 carbons are slightly higher than those with more than 50 branches per 1000 carbons.Differences between the two cooling methods did not affect the value of T g .
On the other hand, our model's T g values fluctuated between 160 and 165 K when the branch content was below 50 per 1000 carbons.The cited literature [53,54,56] also shows a slight decrease in T g with increase in branch content, which is the opposite observation from the model.One of the few factors contributing to this small difference could be the comonomer distribution, controlled by the type of  catalyst used for synthesis or the degree of crystallinity.However, models would have to be systematically constructed to test the effect of comonomer distribution.The LLDPE produced via metallocene catalysis showed a more homogenous comonomer distribution than those obtained through Ziegler-Natta catalysis [53].Generally, T g values of the modelled LLDPE are within an error of the experimentally determined values of commercial polyethylenes.The glass transition temperature of LLDPE was also examined on the same number of branches but at different branch lengths (Fig. 12).It is clear that the length of the branch affects the glass transition temperature for long-chain molecules but not for short-chain molecules.It was important to look at the effect of branch length on glass transition temperature, since published results showed that branch length does affect T g [53,56].Three types of LLDPE (butene-LLDPE, hexene-LLDPE and octene-LLDPE) formed by different α-olefins and ethylene distinguish the branch lengths of this polymer molecule.Unlike  short branches, long branches may influence the crystallisation process and favour the overall crystallinity of LLDPE [56].
Another way to validate whether the models of LLDPE and wax represent the commercial materials was to determine the models' solubility parameters and their behaviour with respect to temperature.The solubility of polymers plays a vital role in creating blends (both melt and solution mixing) and forms part of the molecular interaction theory, which is why the solubility parameter appears in many blend simulation equations.It can also help to simulate conditions suitable for the physical blending of polymers.Therefore, the solubility parameter of LLDPE and wax was determined from the cohesive energy density theory through a module found in Forcite.The solubility of modelled LLDPE made from 50 percent of each monomer decreased observably to a chain length of 180 and then fluctuated at an average value (Fig. 13(a)).The solubility increased from 10 to 30 branches and quickly decreased until 90 branches, where a plateau was reached (Fig. 13(b)).Figure 13(c) clearly shows that the solubility relationship with temperature is almost linear if the best fit is used.Solubility is expected to increase when polymers are heated because the chains become flexible until they are completely mobile.The increase in solubility with increase in temperature agrees with the experimental observations made previously at various temperatures [57].
The solubility parameter has played another special role in atomistic molecular dynamics by forming a bridge between computational cost and the level of accuracy of the models compared to real experimental results.The number of repeat units, chain length and degree of monomer content as a function of solubility have mostly been used to determine the minimum molecular size (of the model), which could represent the real polymer chain length of the blended components [20,58].The reason for finding the minimum molecular size is because the experimentally determined polymer chain length is too big to be handled by the computer simulation software.This method also avoids lengthy dynamic simulations and huge CPU costs [59].Figure 13(b) shows that the solubility parameter remained almost constant for chain lengths of 280 and higher.The solubility of the long chains of LLDPE with different branch contents is presented for three temperatures (Table 2).From 10 to 100 branches, the change in solubility parameter is not consistent with the increase in branch content.However, an increase in temperature shows a clear decrease in solubility, which is a similar trend observed from Choi's models [43,60] and the experimental findings by Han et al. [61].The results in Table 2 complement what has been presented in Fig. 13(c).
The T g values, determined through MD simulations, are from 160 to 180 K.These are relatively close to the experimental values of T g , also referred to as the γ-relaxation, from 147 to 160 K [55].
Models from NPT MD simulations provided information on the changes in volume versus temperature for slowly cooled wax below its melting point (Fig. 14(a)).Most of the constructed wax models showed minor variations in the volume versus temperature plots.Wax is almost 100% crystalline, and it is not expected to show a glass transition.Therefore, the volume change is probably due to volume contraction between the formed wax crystals as the spaces between the chains are reduced during the lamellar packing from the melt.
Figure 14(b) shows the change in density for the modelled wax with different chain lengths plotted against temperature.The same small change in density is seen at 275 K for most wax models.The change in density is also very small, between 275 and 350 K, because the wax chains mostly melt around these temperatures.Density / g cm -3 Temperature / K (b) Because the miscibility of two polymers can depend on the solubility of each component when they are mixed, the effect of temperature on the solubility of the different wax chains and their sizes was investigated (Fig. 15).It was done to establish any relationship between wax solubility and temperature based on the size of the modelled wax chain and the physical state of the wax.There was no apparent trend in the solubility of wax with respect to temperature and the wax chain length between 100 and 275 K.However, the solubility decreased for all the chain lengths when the temperature was risen from 300 to 350 K.At this stage, the van der Waals forces holding the chains were probably weakened by the increase in temperature and the solid wax was converted into a melt.Since there was no observable trend between the solubility parameter and the length of the chains, it did not matter which wax models were used to produce the blend systems with LLDPE.

Conclusions
The method development and validation were successful for both the linear LLDPE and wax models.The modelled wax seemed to have equilibrated well under the chosen MD simulation conditions and parameters.The models' equilibration and validation levels were in line (within a 90% confidence level) with their experimentally determined properties, such as glass transition temperature, solubility parameter, and density, as reported in other published results [53][54][55]62].It should, however, be kept in mind that molecular weight, branch content and length of LLDPE used in experiments are slightly different from the created models.Generally, it was found that (a) the transition temperature of LLDPE was independent of the branch content or the branch lengths, (b) the glass transition temperature of wax was undetectable from the models, which was expected for crystalline materials, and (c) it was difficult to use solubility parameter to find the minimum chain length required to represent a single LLDPE chain molecule's actual size because LLDPE is a copolymer, and its properties are affected by the repeat units in the backbone, as well as the lengths of the branches and their distribution.There was no relationship between the solubility parameter of wax and its chain length.No clear effect of the branch content on the modelled solubility of the LLDPE molecules could be seen.However, the solubility parameter of the LLDPE models decreased when temperature was increased.
The models' properties were validated within a certain degree of error against the polymer' actual physical properties from published data.The different approaches (thermal expansion coefficient, volume and density versus temperature) of determining the glass transition temperature by NPT simulations gave the same results.To a certain level, these developed models represent wax and LLDPE, and they could be used for molecular dynamics simulations, especially for the modelling of the properties, such as solubility and behaviour with respect to temperature, which is important during blend formation.otherwise in a credit line to the material.If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Fig. 1 a
Fig. 1 a Ethene radical molecular structure; b Hexene radical molecular structure (purple colour represent carbon atoms and grey colour represent hydrogen atoms)

Fig. 2 a
Fig. 2 a Short (100 carbons long and 6 branches) and b long (500 carbons long and 70 branches) LLDPE molecular structures (purple colour represent carbon atoms in backbone, dark grey represents carbons in the branches and light grey colour represents hydrogen atoms)

Fig. 3
Fig. 3 Energy change of ethylene and hexene versus the number of optimisation steps

Fig. 4 aFig. 5 a
Fig. 4 a Optimisation energy versus the number of optimisation steps of five different optimisation algorithms; b Zoomed version of (a) showing the optimisation energy between 0 and 1000 steps

Fig. 6 a
Fig. 6 a Energy and pressure of the wax-18 model as a function of MD simulation time; b Density of the wax-18 model as a function of MD simulation time at different temperatures

Fig. 7 Fig. 8
Fig. 7 Simulated density of a LLDPE (long chains / 1000 carbons in the backbone) per different branch content at 430 K; b single LLDPE chains (short chains / 100 carbons in the backbone) based on different branches at 430 K; c wax (at 298 K) in the unit cell (10 chains per unit cell) versus their chain length

Fig. 9 Fig. 10 Fig. 11
Fig. 9 Density as a function of temperature in the modelled LLDPE a with different branch content per 100 carbons in the backbone; b Enlarged section between 120 and 180 K

Fig. 12 a
Fig. 12 a Density versus temperature of the modelled LLDPE with 70 branches per 1000 carbons in the backbone; b Density versus temperature of the modelled LLDPE with 10 branches per 1000 carbons in the backbone

Fig. 13 a
Fig. 13 a Solubility parameter versus the chain length of modelled LLDPE made by a 1:1 ratio of monomers; b Solubility parameter versus the number of branches for the modelled LLDPE (branches per 1000 carbons in the backbone); c Solubility parameter versus tem-

Fig. 14 a
Fig. 14 a Cell volume versus temperature for the wax models with different chain lengths; b Density versus temperature for the wax models with different chain lengths

Fig. 15
Fig. 15 Solubility parameter versus temperature for the modelled wax with different chain lengths

Table 1
Glass transition temperatures calculated for modelled LLDPE by Materials Studio software, and the experimental values from the DMA analysis by Simanke et al. [53]* and Matthews et al. [55]** SC slow cooling, Q quenched cooling, LHM all the LLDPE-1-hexene models in simulations a branch content in mol % b number of branches per 1000 carbons

Table 2
Calculated solubility parameters of the LLDPE at three different temperatures