Young’s modulus, yield point, Poisson’s ratio

Details on the calculation of atomic stresses can be found in the ReaxFF manual .

Note

The calculations are computationally demanding. For optimal performance, a parallel execution on a compute cluster is advised. This can best be done by using the remote job management of the GUI

The simulation may take several days.

Tip

For crystals you can also calculate the Young modulus and Poisson ratio with the elastic tensor.

Overview

This advanced ReaxFF tutorial is based on Radue, Jensen, Gowtham, Klimek-McDonald, King and Odegard, J. Polym. Sci. B, 56, 255-264 (2018). It will demonstrate how to calculate stress-strain curves of an epoxy polymer.

Contents:

  • Setting up a strain rate

  • Calculate the atomic stresses

  • Results: stress-strain curves, Young’s modulus, yield points…

Tip

The molecular structure used in this tutorial was kindly provided by Matthew S. Radue. Other highly cross linked epoxy polymers can effectively be created with a novel bond boost acceleration method (see tutorial).

Setting up

The polymer used in this tutorial is referred to as Tetra(-epoxy) in Ref. 1 since it was build from the tetrafunctional resin epoxy TGDDM (tetra-glycidyl-4,4 0 - diaminodiphenylmethane, Araldite MY 721) and DETDA as hardener. The following cross linking reaction was used to create the polymer structure:

/scm-uploads/doc/Tutorials/_images/stress_strain_epoxy_polymer.png

Start by importing the polymer structure into AMSinput

Click here to download the .xyz file tetra_epoxy.xyz
Import the coordinates in AMSinput:
File → Import Coordinates

Before setting up the strain rate, start with the general molecular dynamics settings. Begin by selecting a force field

Force field dispersion/CHONSSi-lg.ff
/scm-uploads/doc/Tutorials/_images/stress_strain_NPT_settings.png
Click on MoreBtn next to Molecular Dynamics
Select 2000000 as Number of steps (500 ps)
Select 2000 as Sample frequency
/scm-uploads/doc/Tutorials/_images/stress_strain_NPT_settings-2.png

Next we set up the thermo- and barostat

Click on MoreBtn next to Thermostat
Click on the AddButton
Select Berendsen from the menu Thermostat
Set the Temperature to 300.15 K and the Damping constant to 100 fs
/scm-uploads/doc/Tutorials/_images/stress_strain_NPT_settings-3.png
Click on MoreBtn next to MD main options
Click on AddButton next to Barostat
Select Berendsen from the menu Barostat
Set the Pressure to 101000 Pa and the Damping constant to 1500 fs
From the Scale menu, select YZ
/scm-uploads/doc/Tutorials/_images/stress_strain_NPT_settings-4.png

Note

We are applying uni-axial strain, i.e. alongside one lattice vector, and apply the barostat only in the lateral directions. This will allow us to account for Poisson contractions.

Setting up the strain rate

The publication Ref. 1 put a linear strain of 20% over the course of 1 ns. Since we are simulating only 500 ps, for the sake of computational efficiency, the total strain reduces to 10%.

The straightforward way to defining this strain is by defining the strain type to be linear, followed by the final length of the vector at the end of the straining.

Let’s first find and note the length of the a-lattice vector.

Model → Lattice
/scm-uploads/doc/Tutorials/_images/stress_strain_lattice_vecs_length.png

The length of the a-lattice vector a (\(\mid\vec{a}\mid\)) is found to be 37.59793 Å. The endpoint, after stretching it 10%, will then be 41.36 Å. To define this strain, change to the MD deformation panel

Model → MD deformation
Click on AddButton
Enter 2000000 as the final step
Enter 41.36, 0, 0 as the final lengths of the lattice vectors
/scm-uploads/doc/Tutorials/_images/stress_strain_volume_regimes.png

Tip

Putting a 0 as final length tells the program to ignore that particular lattice vector. You can learn about the meaning of each input field by hovering your mouse pointer over it.

Save the input files with a suitable name, e.g. tetra_strain_a, and run the calculation. On 16 cores, it will typically take a bit less than 1 day to compute the trajectory depending on your hardware.

Results

The results can be analyzed either using the AMSMovie GUI module or via a python scripting.

Important

In order to get statistically and physically meaningful results, it’s advised to average the mechanical properties, both over different starting structures as well as all three unaxial strain directions. To obtain quality results as in Ref. 1, an averaging over 5 different polymer structures using 3 strains each (15 trajectories) is required.

Your results may look different from the plots below.

Once the calculation is finished, open the results with AMSMovie:

With the job selected in AMSjobs, click SCM → Movie in the menu bar.

Let’s remove the Energy v.s. time graph, which is not useful for our purposes here:

In AMSMovie, click Graph → Delete Graph

The Young’s modulus can be computed via a linear regression on the Stress-Strain curve. AMSMovie already has a preset for computing the Young’s modulus:

In AMSMovie, select MD Properties → Young modulus → XX
/scm-uploads/doc/Tutorials/_images/stress_strain_amsmovie.png

The graph displays the Stress-Strain data, as well as a smoothed version using a moving average and linear regression lines. The legend presents the coefficients of the linear regression and the corresponding r2 value, enclosed in parentheses.

The Young modulus is the slope of the Stress-Strain linear regression, which in this case is 3.39 * 10^4 bar (or 3.39 GPa).

Note that by this linear regression is performed using the whole strain range (from 0 to 0.1). Following Ref. 1, we will now show how to compute the Young modulus in a smaller strain regime (e.g. only considering the strain range from 0 to 0.03).

Select Graph → Analysis to open the panel showing the analysis options.
Select the Linear regression tab, and in both curves set End to 0.03
Click on OK
/scm-uploads/doc/Tutorials/_images/analysis_details.png

In the plot, you can see that the value of the Young modulus (i.e. the slope coefficient of the linear regression) has changed, since now it’s computing the linear regression only using the values in the 0 to 0.03 strain range. The value of the Young modulus in the 0-0.03 strain range is 5.4 GPa.

In principle, one could calculate the Yield Point using the method outlined in Ref. 1, which involves finding the intersection point between a moving average and the 0.2% offset regression line. However, one would need more statistics / data to get any meaningful result (for example, in Ref. 1, an averaging over 5 different polymer structures using 3 strains each was used)

Tip

In the Analysis Graph options (Menu: Graph → Analysis) you can change the window for the moving average.

Finally, we can compute the Poisson’s ratio (obtained from a linear regression of the transverse strains):

In AMSMovie select MD Properties → Poisson’s ratio → XX.
/scm-uploads/doc/Tutorials/_images/poisson_ratio_amsmovie.png

Since the system is amorphous, the Poisson’s ratio is the average of the two slopes (0.352 in this case).

Literature:

1(1,2,3,4,5,6)

M. S. Radue, B. D. Jensen, S. Gowtham, D. R. Klimek‐McDonald, J. A. King and G. M. Odegard, Comparing the mechanical response of di‐, tri‐, and tetra‐functional resin epoxies with reactive molecular dynamics, Inc. J. Polym. Sci., Part B: Polym. Phys. 56, 255–264 (2018).