M3 Materials modelling: a tension experiment in atomistic detail

Part IB laboratory exercise

Michaelmas Term, 2020, Engineering Laboratory, University of Cambridge, Prof. Gábor Csányi (gc121)

In this exercise you will learn how to model materials on the atomic scale, and use the model to investigate how metals respond under tensile stress. Simulating the response in such detail will allow you an insight into what happens on the smallest scales when a material is under load, how the Poisson ratio arises. At the end of the exrcise, you will see an example of how dislocations are generated and interact, which is the basis for the plastic response of the material. Large scale plasticity is of course more complicated, and depends on the interactions of an immense number of dislocations as well as grain boundaries and more rigid inclusions.

In the course of the exercise, you will:

Create and manipulate atomic scale structural models of crystalline matter
Learn how to compute simple elastic properties of single crystal copper (Cu)
Observe a strain controlled tensile test simulation of a notched bar, and the nucleation of dislocations

This activity introduces you to scientific computing, the use of programming and numerical techniques to model, investigate and learn about a technical subject, in this case materials engineering. There is no extended programming required, and no significant software engineering skills introduced, beyond close reading of documentation.

Working on your own

This is a purely computational laboratory exercise, with no assigned laboratory time. You will be working on your own, but there is a forum on Moodle where you can get guidance and help if you have any technical questions or if you get stuck. It is monitored daily by demonstrators, especially between 2-3pm.

The assigned times in the rota are for marking, and serve as deadlines. You are encouraged to not leave things to the last minute, and you can get marked up for this lab at the DPO helpdesk any time before the deadline assigned to your lab group.

Software Environment

You will be using a python package called Atomic Simulation Environment (ASE), which you will need to download if you want to run this notebook on your own computer. It has excellent tutorials at its documentation page, but the first few exercises below will teach you many of the elements of ASE you will need for this project. Also try to install the matscipy package (which needs a working C++ compiler), if that is successful your simulations will be much faster.

It is recommended that you obtain a Python environment on your own computer such as Anaconda. However, much of this exercise can be completed using a cloud-based notebook server, such as Deepnote, in which case you can try installing the needed packages by using !pip install ase, etc. directly in your notebook. Not all cloud-based servers will let you install packages, and the second part of the exercise is much easier done on your own computer. But Deepnote for example lets you open text terminal (console) to the virtual machine it is running on, so you can execute all the ! commands directly on the console.

The !pip install ase and !pip install matscipy commands also work in an anaconda Python console or a Jupyter notebook that is running on your own computer, or without the "!" in a normal command terminal.

Note: There are several files associated with this notebook, both images and python code (e.g. Morse.py), make sure they are in the same folder as this notebook.

There is also a Morse module inside ASE (under ase.calculators), but that is not implemented correctly (at the time of writing), so make sure below that you import the one that is supplied here and saved in the current directory.

Notes for various operating systems

On Mac OS, it is recommended that you do not work with the Python supplied by Apple, but get your own, e.g. Anaconda, or use Homebrew.
On Windows, you will need MS Visual C++ in order to use the matscipy package (MSVC++ version 14 for use with Python versions > 3.5). Windows is difficult from the point of view of scientific computing, and you may not be able to install matscipy in the end. The lab can be done without it.
If you are working on the Linux workstations in the Department, see this page on how to install new Python packages for yourself.

Software Engineering

Much of this lab is about the application of programming in order to model engineering problem, rather than learning new programming skills. Use the techniques that you learned in the Part IA Lent Computing Exercise, and review them before you start if necessary. In particular:

The first half of the exercise is to be done in this notebook right here.
Please use Python 3 rather than earlier versions.
Keep to the good code documentation practices you learned in Part IA.

It will be assumed from now on that you know how to use these tools.

Materials and structural engineering

Since the topic of this exercise is the response of matter to loading, it will be beneficial if you review the related concepts in the IA Materials and Structures courses, particularly on stress, strain, bulk modulus, shear modulus, Poisson ratio, and dislocations.

Learning objectives and assessment

Learning objectives

Scientific Computing

Making connections between abstract concepts of elastic deformation and numerical representations in a compute model
Ability to calculate elastic properties of materials numerically using an atomistic model
Identifying atomic scale processes underlying plastic deformation via dislocation nucleation and glide

Computer skills

Effective use of python programming, as developed in Part IA
Appropriate use of NumPy arrays, functions, Python classes, class members and member functions.

Assessment guidelines

The following criteria will be used in assessing your implementation. Markers will want to view your Notebook for the elasticity calculations, see your analysis for the tension test simulation.

Your Notebook should display an understanding of the basic aspects of materials engineering involved. Your results should have the correct units, plots have axes correctly labeled.

Code

Your notebook and programs should execute without error.
Programs should be correct and achieve the specified deliverables.
Numerical results should be printed with notes on what they represent, and appropriate units.
Clarity and structure of the code implementations.
Re-use of functions and python objects where appropriate.
Documentation of your functions, and comments in the Notebook and standalone code to enhance readability

Deliverables and Deadlines

The rest of this notebook contains background information that you need to understand the exercise, bits of example code that teaches you how to use the ASE package, and specific deliverables (marked in bold) that you need to present to the markers.

We start with basic usage of ASE, computation of the elastic properties of a single crystal of copper. Do the tasks right here in this notebook. Bring the completed notebook to the a marking session. Do the exercise in your own time, but make sure you get marked by the demonstrators in the DPO by your allocated deadline, which is given for your lab group in the rota. There is a list of marking dates and times here, but these are just the final deadlines, you can go to get marked any time in the DPO when the helpdesk is open, 2-4pm every day.

Task 1.1

You are going to use a simple model to describe the interaction of atoms, called the "Morse potential", which assumes that atoms interact with each other pairwise, and the total potential energy of a collection of atoms is

$$ E{total} = \sum{i < j} V{M}(r{ij}) $$

where the double sum runs over the indices of every pair of atoms (counting every pair only once), $r{ij}$ is the distance between atoms $i$ and $j$, and the _interaction potential is given by

$$ V_{M}(r) = D \lefte^{-2\alpha (r-r_0))} - 2e^{-\alpha (r-r_0))} \right $$

There are three adjustable parameters in this model. You can think of $D$ as fixing the energy scale of the model, $r_0$ as fixing the distance between neighbouring atoms, and $\alpha$ as fixing the length scale over which the interaction between atoms decays (in inverse length units).

In the small scale world of atoms, it is convenient to measure energies in electron Volts (eV), with 1 eV $\approx 1.6 \times 10^{-19}~$J, and distances in Ångstroms (Å), with 1 Å = $10^{-10}~$m, and these are the intrinsic units of ASE. All other derived units follow from this, e.g. forces are in eV/Å, stresses are in eV/Å$^3$. ASE provides the constant GPa that can be used to convert pressures and stresses from eV/Å$^3$ into GPa units.

Here is how to create an object that represents just two copper atoms with a distance 2.5 Å between them:

!pip install ase==3.20.1

from __future__ import print_function import numpy as np from ase import Atoms from ase.units import eV, Ang, GPa d = 2.5*Ang a = Atoms('2Cu', positions=[(0., 0., 0.), (0., 0., d)])

The positions array contains the x, y, z coordinates of each atom, in that order.

Now attach a calculator (ASE-speak for a model that can compute energies, forces and stresses on atoms) to this object, and compute its potential energy. (Any time you create a new atoms object, you need to set its calculator, but you can reuse the same calculator object)

import sys sys.path.append(r'files') import Morse calc = Morse.MorsePotential() a.set_calculator(calc) a.get_potential_energy()

Changing the position of the second atom allows the evaluation of the potential energy as a function of distance

p = a.get_positions() p

p[1,2] = 2.8 p

a.set_positions(p) a.get_potential_energy()

You can also manipulate the positions of the atoms directly by access the array inside the Atoms object:

a.positions[1,2] = 2.9 a.get_positions()

The forces exerted by the atoms on one another can be obtained analogously. Remember that the force is the negative of the 3-dimensional gradient vector of the potential energy. (But there is no directly accessible array in the Atoms object that holds the forces, you have to use the function call)

f = a.get_forces() f

Deliverable 1.1

Write a function that computes the Morse potential energy for two atoms for a given distance between them, and use it to create a plot of the energy against distance. Now do the same for the magnitude of the force exerted by one atom on the other.

What is the distance between the two atoms corresponding to the lowest potential energy ?

Hint: You can create an array of distances (think about what is a sensible range, you don't want atoms to be too close!), and for each one create the Atoms object and compute the potential energy and forces.

Task and Deliverable 1.2

Write a unit test that verifies that the forces returned by the get_forces() function is really the negative gradient of the energy (as returned by the get_potential_energy() function) with respect to the atomic positions. Do this by comparing the returned forces to those obtained by numerical finite differences of energies between two sets of atomic positions, displaced by small amount.

Hint:

Use the idea of the definition of the gradient and its relation to the Taylor expansion:

$$ \nabla f(x) \approx \frac{f(x+\epsilon)-f(x)}{\epsilon} $$

With just two atoms in the "system", you can treat the energy and force as just functions of a scalar variable, the distance between the two atoms. More generally, both the energy and the forces are functions of many atomic coordinates, which we can collect into a vector $\bf{R}$, so the truncated Taylor expansion takes the form:

$$ \nabla f({\bf R})\cdot {\bf d} \approx \frac{f ({\bf R}+\epsilon {\bf d} )-f({\bf R})}{\epsilon} $$

where ${\bf d}$ is the vector representing a small displacement.

Experiment with different values of $\epsilon$ and look at how accurate the approximation is as a function of $\epsilon$. What goes wrong if you make $\epsilon$ too small?

### Importing required packages from __future__ import print_function import numpy as np from ase import Atoms from ase.units import eV, Ang, GPa from matplotlib import pyplot as plt ### Importing and setting up the Morse calculator import sys sys.path.append(r'files') import Morse calc = Morse.MorsePotential() ### Create array of atoms with a range of spacings, as well as delta shifted atoms r0 = 0.5*Ang delta = 1e-11*Ang # for delta < 1e-14, grad_e becomes 0, as the differences are too small for the precision stored in the floats. Also becomes less accurate as you approach this point r_range = np.arange(r0, 10*r0, 0.1*r0) a_range = np.empty(r_range.shape, dtype = object) da_range = np.empty(r_range.shape, dtype = object) for a in range(a_range.shape[0]): a_range[a] = (Atoms('2Cu', positions=[(0., 0., 0.), (0., 0., r_range[a])])) da_range[a] = (Atoms('2Cu', positions=[(0., 0., 0.), (0., 0., r_range[a] + delta)])) ### Calculate potential for each atom pair for a in a_range: a.set_calculator(calc) for a in da_range: a.set_calculator(calc) e = np.empty(r_range.shape, dtype = object) de = np.empty(r_range.shape, dtype = object) for a in range(a_range.shape[0]): e[a] = a_range[a].get_potential_energy() de[a] = da_range[a].get_potential_energy() ### Calculate force between each atom pair f_vec = np.empty(r_range.shape, dtype = object) df_vec = np.empty(r_range.shape, dtype = object) f_mag = np.empty(r_range.shape) df_mag = np.empty(r_range.shape) for a in range(a_range.shape[0]): f_vec[a] = a_range[a].get_forces() df_vec[a] = da_range[a].get_forces() f_mag[a] = abs(f_vec[a][0,2]) df_mag[a] = abs(df_vec[a][0,2]) ### Calculate Derivative grad_e = (de - e) * (1/delta) ### Graph plotting plt.xlabel("distance between centres (Å)") plt.plot(r_range, f_mag, color="tab:red", label = "force on atom (eV/Å)") plt.plot(r_range, -grad_e, ":", color="tab:blue", label = "derivative of Morse potential (eV/Å)") plt.show() ### Similarity test tol = 1e-3 test = True n_test = 0 for i in range(f_mag.shape[0]): x = (f_mag[i]+grad_e[i])/f_mag[i] #print(f_mag[i], -grad_e[i]) ### at some stage, grad_e starts to be the same as f_mag instead of the -ve, not sure why if x > tol: test = False print(test, n_test, f_mag.shape[0], x) break else: n_test+=1 if test == True: print(test, n_test, f_mag.shape[0])

Task 1.3

Create a cubic unit cell of the Cu crystal, and extract some of its properties.

In order to investigate the properties of bulk copper, we need to model a large number of atoms. But evaluating the total potential energy of a large number of atoms takes a long time, and so we need a shortcut. In fact many simple properties of a crystalline solid can be evaluated by just considering its unit cell, i.e. the smallest repeating unit from which the crystal is made. This is typically true for static properties that do not depend on atoms experiencing a variety of neighbour environments, i.e. in the absence of defects. Such properties include the lattice constant (i.e. the density), the equation of state, the elastic constants (including the bulk modulus), the Poisson ratio, etc.

The trick is to assume that the atoms are arranged in perfect crystalline order, and only explicitly consider those atoms that are in a single unit cell. In order to correctly evaluate the energy, the effect of atoms in neighbouring unit cells needs to be included, but this can be done by using periodic boundary conditions. The page behind the link is rather detailed, but think about the Fourier Series from IA maths, which was applicable to periodic functions, but you only ever had to consider the function values over a single period. This is similar.

In the image below, the blue shaded area is the unit cell, its sides are the lattice vectors that generate the periodic images.