Artificial Neural Networks for Predicting Mechanical Properties of
  Crystalline Polyamide12 via Molecular Dynamics Simulations
Predicting material properties of 3D printed polymer products is a challenge
in additive manufacturing due to the highly localized and complex manufacturing
process. The microstructure of such products is fundamentally different from
the ones obtained by using conventional manufacturing methods, which makes the
task even more difficult. As a first step of a systematic multiscale approach,
in this work, we have developed an artificial neural network (ANN) to predict
the mechanical properties of the crystalline form of Polyamide12 (PA12) based
on data collected from molecular dynamics (MD) simulations. Using the machine
learning approach, we are able to predict the stress-strain relations of PA12
once the macroscale deformation gradient is provided as an input to the ANN. We
have shown that this is an efficient and accurate approach, which can provide a
three-dimensional molecular-level anisotropic stress-strain relation of PA12
for any macroscale mechanics model, such as finite element modeling at
arbitrary quadrature points. This work lays the foundation for a multiscale
finite element method for simulating semicrystalline polymers, which will be
published as a separate study.
Caglar Tamur, Shaofan Li, Danielle Zeng
Artificial neural networks for predicting mechanical properties of crystalline polyamide12 via molecular dynamics simulations

###### Abstract

Predicting material properties of 3D printed polymer products is a challenge in additive manufacturing due to the highly localized and complex manufacturing process. The microstructure of such products is fundamentally different from the ones obtained by using conventional manufacturing methods, which makes the task even more difficult. As a first step of a systematic multiscale approach, in this work, we have developed an artificial neural network (ANN) to predict the mechanical properties of the crystalline form of Polyamide12 (PA12) based on data collected from molecular dynamics (MD) simulations. Using the machine learning approach, we are able to predict the stress-strain relations of PA12 once the macroscale deformation gradient is provided as an input to the ANN. We have shown that this is an efficient and accurate approach, which can provide a three-dimensional molecular-level anisotropic stress-strain relation of PA12 for any macroscale mechanics model, such as finite element modeling at arbitrary quadrature points. This work lays the foundation of a multiscale finite element method for simulating semicrystalline polymers, which will be published as a seperate study.

keywords: _Molecular dynamics, crystalline polymers, rexaff, polyamide, machine learning, neural network_

## 1 Introduction

Material modeling of 3D printed material products has been a great challenge for additive manufacturing and other advanced manufacturing technologies. In additive manufacturing, the product is built on a material point upon a material point; thus, the resultant properties are highly dependent on the printing process both spatially and temporally at the localized regions. A microstructurally informed and robust mechanical model is highly anticipated, as it would help to overcome the challenges with manufacturing process, product design and material selection.

One of the promising approaches is the multiscale modeling of 3D printed materials [1; 2; 3], which has the ability to extract bottom-up material information to assess the
macroscale properties of a product that is additively manufactured. However, the cost of multiscale modeling may vary, especially when it is involved with molecular-scale simulations. Traditionally, in multiscale modeling of crystalline materials, this is overcome by using a technique called the Cauchy-Born rule [4, 5], which is an approximation of molecular dynamics by using simplified molecular statics. However, for semicrystalline or amorphous materials, due to the lack of a definite crystal structure, full-scale molecular dynamics must be utilized, which will significantly increase the cost of the multiscale simulation. This is because the macroscale finite element method would require a molecular level stress-strain relation at each quadrature point for an arbitrarily given strain, and it is almost impossible to conduct concurrent molecular dynamics simulations unless one has abundant computational resources.

The task at hand is to simulate the molecular-level stress-strain relation under arbitrary strain and temperature with good accuracy, efficiency, and low cost. To solve this problem, in the present work, we adopted a machine learning approach by developing an artificial neural network (ANN) that is trained on a molecular dynamics (MD) dataset, which can predict the molecular-level mechanical material properties of polymeric materials.

Current efforts in machine-learning driven material property predictions involve a variety of approaches. One common approach is to use the features based on the material composition to obtain structure and electronic properties, such as predicting bandgap of semiconductor materials using convolutional neural networks [6] and predicting adsorption energies of high-entropy alloys through deep neural networks [7]. In a recent study, standard feed-forward neural networks are used with EAM based MD simulations to model crystal plasticity in a multiscale framework [5].

Polymeric materials have three different microstructural forms: crystal, amorphous, and semicrystalline. As a starting point, we shall first develop an artificial neural network to predict the mechanical response of the anisotropic crystalline forms of polymers. The approach we presented can be extended to include amorphous phases and would eventually incorporated into the constitutive relations of a multiscale finite element method for semicrystalline polymers.

One of the polymeric materials widely used in additive manufacturing is Polyamide12 (PA12), also known as Nylon12, which is a synthetic polymer that is used in many industrial applications such as automobile parts, aerospace applications and medical components. An automotive part manufactured using the Multi Jet Fusion process is shown in Fig. 1.
PA12 has the chemical formula [-(CH\({}_{2}\))\({}_{11}\)C(O)NH-]\({}_{n}\) as visualized in Fig. 2. In recent years, many 3D printed PA12 products have been fabricated, and they have shown some outstanding material properties. Additively manufactured PA12 is a semicrystalline material, in which the microstructure has a sandwich-like structure of alternating regions of crystal and amorphous zones.

Its structural characterization has been extensively studied since the 1970s with x-ray diffraction (XRD) experiments [8; 9] and with nuclear magnetic resonance (NMR) spectroscopy [10], which revealed that PA12 displays polymorphism and can have crystal phases \(\alpha,\alpha^{\prime},\gamma\) and \(\gamma^{\prime}\). The most abundant phase is the \(\gamma\) form, which results from slow cooling at atmospheric pressure, whereas the other phases require specific conditions such as rapid quenching and/or high-pressure treatment. For the purposes of our study, we focus on the \(\gamma\) crystal form.

Figure 1: 3D printed Nylon PA12 auto part via HP multi-jet fusion (Courtesy of Ford Motor Company).

Figure 2: Polyamide12 molecules. Gray, white, blue and red spheres represent carbon, hydrogen, nitrogen and oxygen atoms respectively. Dashed lines indicate the H bonds.
## 2 Molecular dyannics simulations

In this section, we describe the details of molecular dynamics simulations, such as the preparation of initial configurations, the force field selection, and simulation procedures, and we conduct a preliminary study.

### System setup

We start by constructing the unit cell of the \(\gamma\) form PA12 crystal. The \(\gamma\) phase has a pseudo-hexagonal monoclinic structure with the lattice parameters summarized in Table 1[8; 9]. Using the atomistic coordinates and the lattice parameters from the experimental literature [9], the unit cell of the \(\gamma\) phase is constructed. This monoclinic structure contains four PA12 chains and is visualized in Fig. 3. Note that the dashed lines indicate hydrogen bonds, and the unit cell is periodic in all directions.
### Force field selection

Force fields represent the interactions between atoms and molecules using a set of equations. In classical forms, force fields consist of empirical potentials that describe bonded interactions (primary bonds, bond angles, and dihedrals) and non-bonded interactions (van der Waals and electrostatic). Choosing a suitable force field for the system under consideration is a crucial part of molecular dynamics simulations.

Several empirical force fields have been used in the simulation of polyamides; some popular examples being OPLS, CVFF, COMPASS, and DREIDING [11]. However, these potentials are typically unable to model chemical reactions and cannot account for bond dissociation. Hence, they are not the ideal choice for investigating the mechanical response of the system in the context of fracture mechanics and rupture of polymer crystals, which may involve the chain scission process [12; 13].

Reactive force fields, on the other hand, have been developed to simulate complex chemical reactions and have been used successfully in deformation simulations that involve bond cleavage. Once such a force field is ReaxFF, originally developed by van Duin [14], which utilizes the so-called bond order formalism to describe chemical bonding. ReaxFF requires no topology information, and it is easy to construct initial MD configurations since only the atomistic coordinates are needed. ReaxFF specially treats weak hydrogen bonds with an explicit energy term in the functional, which plays an important role in polymer systems. The current implementation in LAMMPS utilizes the functional described in [15].

ReaxFF is parameterized for a wide range of materials through quantum mechanical calculations [16] and has also been used in the simulation of polyamide crystals before [17]. For our polymer system, we tested three such parameterizations [18; 19; 20]. The test process involved relaxing the initial structure in the NPT ensemble and comparing the resultant lattice parameters with the experimental values. Among these, the ReaxFF parameterization by Mattson [19] yielded the best result; therefore, it was chosen as the force field for this study.

Figure 4: Supercell of perfect PA12 crystals (a) View from x-y plane (b) View from y-z plane.
### Simulation details and preliminary analysis

Molecular dynamics simulations are performed using the Large-Scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) [21]. The simulations are run on Berkeley's high-performance computing cluster Savio using a \([4\times 2\times 4]\) CPU grid. MD time steps are chosen as 0.5 fs. Nose-Hoover thermostat and barostat [22] are used to control temperature and pressure, respectively, with damping parameters chosen as 50 fs and 500 fs.

The initial structure is relaxed in the constant pressure and temperature (NPT) ensemble at 300 K and 1 atm, where the pressure is controlled independently in all directions. We observed that an NPT simulation for 50 ps was sufficient to reach equilibrium, as the lattice parameters stabilized adequately at this point.

As a preliminary study and to investigate in detail the anisotropic mechanical response of PA12 crystals, we constructed a larger MD cell with 38,000 atoms. After relaxation under ambient conditions, we deform the system in different directions under uniaxial tension. Uniaxial deformation is imposed by stretching the unit cell in one direction at a constant strain rate and relaxing the other two dimensions with the NPT ensemble at 1 atm. For this part of the study, we used slower strain rates and deformed the cell quasi-statically in a step-by-step fashion as described in [23]. In this process, the MD cell is deformed to 1% of its final stretch, relaxed in the NVT ensemble, and the stress is sampled and averaged over 5 ps intervals. These steps are repeated until rupture of the polymer chains is observed and the resultant stress-strain behaviors are shown in Fig. 5 for each direction.

From Fig. 5, we observe that the system is highly anisotropic; crystal is significantly more stiff in the z-direction, which is the direction of polymer chain alignment, and it is more ductile in the other two directions. Interestingly, for the x-direction Fig. 5 (a), there was a remarkable increase in the elastic modulus once the strain reaches around 12%. We explored this phenomenon to comprehend whether it involves a mechanically induced phase transformation, by deforming the original monoclinic crystal supercell in the x-direction. Resulting atomistic configurations before and after the change in elastic modulus are visualized in Fig. 6, in which we clearly observe a change in the microstructural configurations.

Figure 5: Stress (MPa) versus Strain (engineering) behavior of PA12 in uniaxial tension.
Upon additional analysis, we understood that the process is completely reversible and the effect disappears once the MD cell is unloaded to its original state. Further investigation may be required to draw a meaningful conclusion on the phenomenon.

Figure 6: Snapshots of the atomistic configurations during deformation in the x-direction.
## 3 Artificial neural network for constitutive law

In this section, we present a data-driven approach to model the hyperelastic constitutive law of crystalline PA12, to be used in multiscale mechanics simulations. Our task is to find the mapping between the deformation state and the resulting material response, using the data obtained from the molecular dynamics simulations on an artificial neural network (ANN). We discuss data collection, model selection, training phase, and predictions in detail.

### Data collection and processing

In order to generate the data set to train the learning model, we performed MD simulations by deforming the PA12 supercell in different directions with varying strain rates. For the purposes of this study, we limit our attention to uniaxial and biaxial tensile deformations, although the presented methodology should be applicable to any deformation mode. Using the findings of the preliminary study, Fig.5, we determined the final stretch limits of the simulation box necessary in each direction. Accordingly, we chose the ranges of strain rates in each direction, which are summarized in Table 2. In all MD simulations, the supercell is deformed for 10,000 steps with a rate chosen from the prescribed ranges depending on the loading direction. Specifically, we chose 15 evenly spaced values from the ranges in Table 2 and took their combinations to obtain uniaxial and biaxial loading cases, resulting in 720 different MD simulations each having different final stretches and strain rates.

Each deformation simulation is run for 5 ps and the Virial stress and the box dimensions are sampled at 50 fs intervals. Consequently, after running all the 720 simulations, we collected a data set consisting of 72,720 data points. Each data point is 12 dimensional, which represents the deformation (6 components) and stress (6 components) states of the PA12 crystal. Note that the process described here involves remarkably high strain rates and short simulation times, since it was not computationally feasible to conduct slower simulations and collect a large data set at the same time.

The original data from MD simulations consist of box dimensions \(l=\{l_{x},l_{y},l_{z}\}\) which can be used to construct the deformation gradient \(\mathbf{F}\) and the stress measure we get is the pressure tensor, known as the Virial stress. Virial stress has been shown to be equivalent to continuum Cauchy stress \(\mathbf{\sigma}\)[24]. In order to adopt the model into the continuum mechanics framework, we need to transform the data set into the tensor quantities utilized in finite deformation. To preserve material objectivity [25; 5], we chose the energetic conjugates right Cauchy-Green tensor \(\mathbf{C}\) and the second Piola-Kirchhoff (PK2) stress tensor \(\mathbf{S}\) to represent the deformation state and material response, respectively. Each of these second-order tensors have nine components, but making use of the symmetry and Voigt notation we represent

\begin{table}
\begin{tabular}{c c c c} \hline \hline  & \multicolumn{3}{c}{Deformation Direction} \\  & x & y & z \\ \hline Ultimate Strain & 0.5 & 0.6 & 0.15 \\ Strain Rate (\(10^{-6}/fs\)) & [3.3, 50] & [4; 60] & [1; 15] \\ \hline \hline \end{tabular}
\end{table}
Table 2: Ultimate strain and resulting ranges of strain rates for each direction.
them in the vector form as below.

\[\mathbf{C} =[C_{11},\ C_{22},\ C_{33},\ C_{12},\ C_{13},\ C_{23}] \tag{1}\] \[\mathbf{S} =[S_{11},\ S_{22},\ S_{33},\ S_{12},\ S_{13},\ S_{23}] \tag{2}\]

Using the elementary relations of continuum mechanics, we can obtain the \(\mathbf{C}\) and \(\mathbf{S}\) tensors as follows.

\[J =\det\mathbf{F} \tag{3}\] \[\mathbf{C} =\mathbf{F}^{T}\mathbf{F}\] (4) \[\mathbf{S} =J\mathbf{F}^{-1}\mathbf{\sigma}\mathbf{F}^{-T} \tag{5}\]

Now we can state our goal as finding the map,

\[\Psi:\mathbf{C}\rightarrow\mathbf{S} \tag{6}\]

where \(\Psi\) encapsulates the constitutive model that we are going to approximate through supervised learning.

Finally, we normalize our input data set to have a mean of zero and a standard deviation of one, a process known _standardization_. It is common practice in gradient-based learning methods to standardize the data, and it improves the performance of ANNs by helping to solve the convergence issues of the backpropagation algorithm [26; 27]. Standardization can be described as follows.

\[C_{ij}^{std}=\frac{C_{ij}-\hat{\mu}_{ij}}{\hat{\sigma}_{ij}}\quad\text{for} \quad\mathrm{i,j}=1,2,3. \tag{7}\]

where \(C_{ij}^{std}\) are the normalized components of the right Cauchy-Green tensor, \(\hat{\mu}_{ij}\) is the sample mean and \(\hat{\sigma}_{ij}\) is the sample standard deviation of the respective component.

### Model selection and results

We adopt a feedforward, fully connected network architecture to construct our regression ANN, which may be called a multilayer perceptron (MLP). As shown by the universal approximation theorem, feedforward neural networks can approximate any continuous function, provided that they have at least one hidden layer, have enough neurons, and the activation functions satisfy certain properties [28; 29]. Thus we believe that a MLP is a suitable choice to approximate the hyperelastic constitutive law. Keras [30] Python interface of the TensorFlow [31] machine learning package is used to select and train our ANN model.

Our task is to find an approximation to the constitutive model defined in Eq.(6). Formally, we can express the learning problem as follows.

Find

\[\mathcal{N}(\mathbf{C},\mathbf{w})=\mathbf{\hat{S}} \tag{8}\]
such that

\[\mathbf{w}=\operatorname*{argmin}_{\mathbf{w}}\mathcal{L}(\mathbf{S},\mathbf{\hat{S}}) \tag{9}\]

where \(\mathcal{N}\) is encodes the ANN, \(\mathbf{w}\) are known as the _weights_ of each neuron and \(\mathcal{L}\) is a _loss function_. The task in supervised learning is to find optimal weights \(\mathbf{w}\) such that the metric defined in the loss function is minimized.

To select the ANN model parameters known as _hyperparameters_, such as the number of hidden layers, number of neurons at each layers, type of activation functions and learning rate of the gradient descent, we conduct hyperparameter optimization. We chose two candidate activation functions which are commonly used in regression, _ReLu_ function and its smooth approximate version known as the _Softplus_ function.

\[\text{ReLu}\ :\ \sigma(\text{x})=\text{max}(\text{x})\qquad\quad\text{ Softplus}\ :\ \sigma(\text{x})=\log(1+\text{e}^{\text{x}}) \tag{10}\]

For the loss function, we choose the mean squared error (MSE) of the PK2 stress, which is defined as

\[\mathcal{L}_{\text{MSE}}(\mathbf{S},\mathbf{\hat{S}})=\frac{1}{n}\sum_{i=1}^{n}||\mathbf{ S}_{i}-\mathbf{\hat{S}}_{i}||_{2}\, \tag{11}\]

where \(\mathbf{S}\) is the PK2 stress obtain from MD simulations, \(\mathbf{\hat{S}}\) is the predicted stress from ANN and \(n\) is the number of data points.

We tried two methodologies to select the best set of hyperparameters, the HyperBand algorithm [32] and Bayesian Optimization with Gaussian Process [33]. Both algorithms would pick the best parameters from a predefined set, which would give the lowest validation loss. The set of parameters that we make the search is determined by the preliminary analysis, where the optimizer of the ANN is chosen as the Adam algorithm [34], which implements a version of the stochastic gradient descent (SGD) method. Accordingly, we chose the following range of parameters to perform hyperparameter tuning.

We start by partitioning the data into a random 80% - 20% train-test split. Then, we run hyperparameter optimization using Hyperband and Bayesian Optimization methods in the domain described in Table 3, leaving aside 20% of the training data to be used to compute validation loss. For all models, the learning rate \(\gamma=10^{-2}\) and the number of hidden layers of four or five gave us the best results. For the rest of the hyperparameters, we investigated the top ten models and chose the ones with the lowest complexity to reduce overfitting. Resulting candidate models are summarized in Table 4.

\begin{table}
\begin{tabular}{c c c c} \hline \hline  & \multicolumn{3}{c}{Hyperparameters} \\ \hline Hidden Layers & Neurons & Activation Function & Learning Rate \\ \([2,5]\) & \([32,256]\) & \(\{\text{ReLu, Softplus}\}\) & \(\{10^{-1},10^{-2},10^{-3}\}\) \\ \hline \hline \end{tabular}
\end{table}
Table 3: Hyperparameter search grid
The models summarized above are trained on the train data set for 1000 epochs, and predictions are made on the test data to assess our final performance. The resulting train-test history curves are presented in Fig. 7. We conclude that the ANNs are stable, there is no significant overfitting, and they perform well against the test data set, as seen from the test loss. The best model proves to be ANN-4, the rightmost column in Table 4, and its architecture is schematically visualized in Fig. 8.

\begin{table}
\begin{tabular}{c c c c c} \hline  & Hyperband (1) & Bayesian (2) & Hyperband (3) & Bayesian (4) \\ \hline Input Layer & 6 & 6 & 6 & 6 \\ Hidden Layer 1 & 192 \(\times\) ReLu & 64 \(\times\) ReLu & 160 \(\times\) ReLu & 224 \(\times\) ReLu \\ Hidden Layer 2 & 32 \(\times\) Softplus & 32 \(\times\) Softplus & 64 \(\times\) Softplus & 96 \(\times\) Softplus \\ Hidden Layer 3 & 32 \(\times\) ReLu & 256 \(\times\) Softplus & 32 \(\times\) ReLu & 64 \(\times\) ReLu \\ Hidden Layer 4 & 64 \(\times\) ReLu & 32 \(\times\) ReLu & 128 Softplus & 96 \(\times\) Softplus \\ Hidden Layer 5 & 128 \(\times\) ReLu & 64 \(\times\) ReLu & - & - \\ Output Layer & 6 \(\times\) Linear & 6 \(\times\)Linear & 6 \(\times\) Linear & 6 \(\times\)Linear \\ Validation Loss & 15630 & 15026 & 14109 & 13538 \\ \hline \end{tabular}
\end{table}
Table 4: ANN architectures resulting from Hyperband method and Bayesian Optimization

Figure 7: Learning curves of the ANNS: Training and test loss during training phase.
Finally, we tested ANN-4 against a new set of uniaxial tension MD simulations to observe our models performance against the variance in strain rate. Two of the simulations utilized strain rates \(\dot{\epsilon}_{1}=5\times 10^{-6}/fs\) and \(\dot{\epsilon}_{2}=10\times 10^{-6}/fs\), that are within the training interval defined in Table 7. Predictions of ANN-4 is compared with the MD results and the resultant \(C_{33}\) versus \(S_{33}\) (MPa) relations are shown in Fig. 9. Predictions in the elastic range and the ultimate strength display great agreement with the MD findings, which proves the predictive power of the ML approach.

Figure 8: ANN-4 schematic representation [35].

Figure 9: Predictions of ANN-4 for uniaxial tension in z-direction for a new set of simulations. The strain rates are contained in the training range.
We consider two additional MD simulations that employed strain rates lying beyond the training range, encompassing both extremes. These strain rates are defined as \(\dot{\epsilon}_{3}=0.5\times 10^{-6}/fs\) (slower) and \(\dot{\epsilon}_{4}=30\times 10^{-6}/fs\) (faster). The corresponding outcomes are depicted in Fig. 10, alongside the MD results, for comparison. While the slower rate led to a slightly higher ultimate strength, a remarkable agreement is observed within the elastic range. Here, the ANN effectively captured intricate aspects of the loading path. The findings clearly demonstrate the model's capacity for generalization and robustness across various strain rates, which makes it a promising candidate for effectively describing the constitutive laws of crystalline polymers to be used in large deformation FEM simulations.

## 4 Conclusion and future perspectives

In this study, we have developed an artificial neural network to model the mechanical response of crystalline polymers based on molecular dynamics simulations. Initial configurations are generated to model PA12 in \(\gamma\) crystal form. ReaxFF was our choice for the force field, and MD simulations are performed under uniaxial and biaxial tensile deformations to generate a data set for the learning model.

ANN architecture is selected through hyperparameter tuning algorithms and carefully calibrated using the training set. The training and validation process produced the best model to be ANN-4, which is proven to be robust and accurate with no significant overfitting. The best model is tested against different strain rates, which covers above and beyond the loading rate range of the training set. The results demonstrated the predictive power and generalization capacity of our approach, even with strain rates that are much slower or faster than the original data.

Figure 10: Predictions of ANN-4 for uniaxial tension in z-direction for a new set of simulations. The strain rates are beyond the range of the training data.
The presented methodology proved to be an efficient way of modeling the microstructurally informed constitutive relation of crystalline polymers, which would be used in macroscale FEM simulations. Namely, we presented an approach to accurately predict the mapping between right Cauchy-Green (\(\mathbf{C}\)) and second PK2 stress (\(\mathbf{S}\)) tensors, that shall be utilized to compute the strain energy density of large deformation problems. Once the model is trained, we can simply treat it as a black box for computing the strain energy at arbitrary quadrature points and time steps, corresponding to the current deformation state of the material.

There are some caveats that need further attention. As detailed in Section 2.3, there may be a mechanically induced phase transformation of PA12, which may depend on the particular additive manufacturing method and process parameters such as temperature and printing speed. Furthermore, we only considered uniaxial and biaxial tensile deformations in our training set. For general loading conditions, such as pure shear, compression, and torsion, the artificial neural network will need more data and the associated ML training procedure may become more complicated and expensive.

It is worth noting that our current investigation was limited to the fully crystalline form of PA12 for purposes of demonstration. The actual material components consist of semicrystalline structures wherein the microstructure comprises both amorphous and crystalline regions. However, the presented MD and ML methodologies can easily be adopted to amorphous domains. Additionally, considering larger simulation cells, having longer relaxation windows and adopting slower strain rates may be necessary to simulate more accurate and realistic scenarios. These matters will be subjected further investigated and will be reported in a separate study.

## 5 Acknowledgements

We acknowledge the support from the University Grant by Ford Motor Company. Caglar Tamur gratefully acknowledges the financial support for his graduate studies and research by the Fulbright Scholarship Program.