Part II: Workflows in Materials Science #

Now that we are familiar with python, we will look a typical workflow in materials science. We will calculate the Youngs modulus, first with LAMMPS, a widely used simulation code. If you are looking for a quick way to calculate elastic constants, example script provided with LAMMPS is the way to go. Here, the method itself will be explored. This method will be followed here and is a good resource for further reading.

A crucial input that we need is a way to model to the interaction between the atoms. Here will we will use an Embedded Atom Method interatomic potential to simulate Aluminium.

Using Stress#

For a cubic material, there are three independent elastic constants - \(C_{11}\), \(C_{12}\) and \(C_{44}\). \(C_{11}\) and \(C_{12}\) can be calculated from the stress-strain relationships. If a strain \(\epsilon_{11}\) is applied and the stresses \(\sigma_{11}\) and \(\sigma_{22}\) are measured, the elastic constants can be calculated directly from,

\[\begin{split} C_{11} = \frac{\sigma_{11}}{\epsilon_{11}} \\ C_{12} = \frac{\sigma_{22}}{\epsilon_{11}} \end{split}\]

Here, we will focus on \(C_{11}\), which is the Youngs Modulus.

The LAMMPS script is already prepared for you. Let us take a look first.

! cat lammps.in

Now we can run the script

! lmp_serial -in lammps.in

We can read in the output we need, from stress.dat

import numpy as np
import matplotlib.pyplot as plt

lx, sigma11, sigma22 = np.loadtxt("stress.dat", unpack=True)

The first thing to do to look at our data is to convert the simulation box dimension along x to strain. The value at position 0, lx[0] is the undeformed box dimension. Thus, the strain is given as,

strain = (lx -lx[0])/lx

Great! Now we can see the stress values that we calculated,

sigma11

The first thing to notice is that the values are negative. Since we are calculating the pressure, the negative sign just indicates that the simulation box wants to contract. The values are in bar, so to convert to GPa, it is multiplied with 0.0001

sigma11 = -1*0.0001*sigma11
sigma22 = -1*0.0001*sigma22

Now we can plot the strain versus the stress.

plt.plot(strain, sigma11, '-', label=r"$\sigma_{11}$")
plt.plot(strain, sigma22, '-', label=r"$\sigma_{22}$")
plt.xlabel(r"$\epsilon$");
plt.ylabel(r"$\sigma$ GPa");

The slope of \(\sigma_{11}-\epsilon_{11}\) and \(\sigma_{22}-\epsilon_{11}\) will give \(C_{11}\) and \(C_{12}\) respectively. To calculate the slope, we can fit a straight line using numpy. We can use the polyfit function for that. polyfit stands for polynomial fit, and we need to provide an order of fit.

np.polyfit(strain, sigma11, 1)

The function above fits it to \(y = m*x + c\), the output is the slope, \(m\) and the intercept \(c\). Therefore the Youngs modulus is 105 GPa. The expected value would be 118 GPa, so our estimate is fairly close.

Task

Calculate C12

We successfully managed to calculate the elastic constants. However, there were a number of things we had to do or know:

  • How to write LAMMPS script

  • Interatomic potential file

Additionally, we ran LAMMPS from the terminal, then took the results file, and post-processed it with python. These workflows come with a number of risks which affect the reproducibility. Also, how can we be sure that we ourselves can understand and interpret these calculations as time goes by.

  • We need to organise folders and files ourselves (possibly with a naming convention)

  • We need to be able to look into the LAMMPS file and understand the input quantities.

  • We need to save the python script for the analysis we did.

  • ..

This is where workflow environments come in. Workflow environments automate a number of these tasks and make it easy for you to focus on the science. One such workflow environment in pyiron. Now we will do the same calculation using pyiron.

from pyiron import Project

The Project object introduced below is central in pyiron. You can imagine it to be equivalent to the current scientific project you are working on. In this case, we will call it Al_elastic_constant

project = Project("Al_elastic_constant")

The next major object in pyiron is Job. A Job is a class designed to do a particular task at hand. pyiron comes with a lot of pre-built jobs. However, it is also easy to create your own.

In this case we will use a job for calculating elastic constants. You can always press Tab key for autocompletion to see what options are possible. Now we create a LAMMPS job.

job = project.create.job.Lammps("job_1")

First thing we need is a structure! pyiron also provides tools to create and visualise structures.

structure = project.create.structure.bulk('Al', cubic=True)

structure.plot3d()

We assign the structure to a job

job.structure = structure

Next thing we need is an interatomic potential. We can check which potentials are available.

job.list_potentials()

There are a large number of potentials available for Al. We will select the same one as before. However, here we did not have to find the file, or consider how it is mapped into the LAMMPS input file.

job.potential = "2008--Mendelev-M-I--Al--LAMMPS--ipr1"

As we discussed before, the first step is a relaxation or minimization. We will do the same here.

job.calc_minimize()

Now we use the LAMMPS job as a reference, and now tell pyiron to calculate the elastic constants.

elastic_job = job.create_job(project.job_type.ElasticMatrixJob, "elastic_job")

We provide the input strain range, 0.001

elastic_job.input["eps_range"] = 0.001

And finally call run

elastic_job.run()

The elastic_job starts several other jobs automatically. Once it is finished, we can look at the results.

elastic_job["output/elasticmatrix"]["C"][0,0]

You can see that we got very similar results. However, we could directly focus on the scientific aspects of the workflow. pyiron takes care of managing the folders, giving IDs to your jobs, and the post-processing automatically

Further reading..