Introduction to Atomistic Simulations with Pyiron

Before the excercise, you should:

• Be familiar with python especially with numerical libraries like numpy and plotting tools like matplotlib

• Understand how Jupyter Notebooks work

The aim of this exercise is to make you familiar with:

• A general overview of what pyiron can do

• How to set up atomic structures and run atomistic simulation codes through pyiron

Jupyter Crash Course

1. Select cells by clicking on them.

2. Navigate through with up and down keys (or k and j for you vimmers).

3. Press Enter to edit a cell.

4. Press Shift-Enter to execute it.

5. Create new cells above or below the current one with a or b.

6. Copy, cut and paste them with c, x and v.

7. Press m to turn a cell into a markdown cell.

8. See the Help in the toolbar for more.

For Today

Importing necessary libraries

As a first step we import the libraries numpy for data analysis and matplotlib for visualization.

from pyiron_atomistics import Project
import matplotlib.pyplot as plt
import numpy as np

Fundamentally, we only need to import one module from pyiron: the Project class

The Project object introduced below is central in pyiron. It allows to name the project as well as to derive all other objects such as structures, jobs etc. without having to import them. Thus, by code completion Tab the respective commands can be found easily.

We now create a pyiron Project named ‘first_steps’.

Creation of a project instance

pr = Project("first_steps")

The project name also applies for the directory that is created for the project. All data generated by this Project object resides in this directory.

pr.path

pr

The groups and nodes will be populated later, as we add jobs and sub project to it.

Creating atomic structures

Every atomistic simulation needs an atomic structure. For more details on generating and manipulating structures, please have a look at our structures example. In this section however, we show how to generate and manipulate bulk crystals, surfaces, etc. pyiron’s structure class is derived from the popular ase package and any ase function to manipulate structures can also be applied here.

Creating a bulk fcc cubic unitcell

Al_unitcell_cubic = pr.create.structure.bulk('Al', cubic=True, a=4.04)
Al_unitcell_cubic

Creating a super cell.

Al_supercell_3_3_3 = Al_unitcell_cubic.repeat([3, 3, 3])
Al_supercell_3_3_3.plot3d()

Running an atomistic calculation using interatomic potentials (with LAMMPS)

Once we have an atomic structure, we can set up a simulation “job” of any atomistic simulation that is intergrated within pyiron. In this section, we are going to use the popular LAMMPS code.

# Create a job
job_lammps = pr.create.job.Lammps(job_name="lammps_job")

Every atomistic simulation code needs an input atomic structure. We use the Al supercell structure we created earlier

# Assign an atomic structure to the job
job_lammps.structure = Al_supercell_3_3_3

Once the structure is assigned, an appropriate potential should also be chosen. This list of available for the structure containing Al can be found below. This list originates from the NIST Interatomic Potential Database.

# See available potentials
job_lammps.list_potentials()[20:30]

# Choose one of these potentials
job_lammps.potential = "2005--Mendelev-M-I--Al-Fe--LAMMPS--ipr1"

At this stage, the computational parameters for the simulation needs to be specified. pyiron parses generic computational parameters into code specific parameters allowing for an easy transition between simulation codes

# specify calculation details: in this case: MD at 800 K in the NPT ensemble (pressure=0) for 10000 steps
job_lammps.calc_md(temperature=800, pressure=0, n_ionic_steps=10000)

We can now see how pyiron sets-up the corresponding LAMMPS input

job_lammps.input.control

Once the run() commmand is called, pyiron creates necessary input files, calls the simulation code, and finally parses and stores the output.

job_lammps.run()

When printing the project, the saved job will also appear under nodes now.

pr

You can get a quick overview with the job_table method.

pr.job_table()

Analysing a calculation

pr

Indexing into a Project loads a job from storage.

job_loaded = pr['lammps_job']

Jobs in turn also have groups and nodes in them, corresponding to the input and output that was stored for them.

job_loaded

You will generally find all relevant data in the output group, for your further processing.

job_loaded["output"]

The group structure holds the final structure at the end of the simulation, whereas the generic group keeps output you can expect for all atomistic jobs in pyiron.

job_loaded["output/generic"]

plt.plot(job_loaded["output/generic/energy_pot"])

Structures from the MD trajectory can be accessed using get_structure, passing as the argument the index of the step you want to retrieve.

final_struct = job_loaded.get_structure(frame=5)
final_struct.plot3d()

Using the plot3d method you can visualize the structure interactively. To see the full trajectory, you can use the animate_structure method.

job_loaded.animate_structure()

Let’s go through a few more examples of visualizing the output of the job.

When you run the job, you can specify how often pyiron and lammps should be writing the intermediate structures and energies, etc. The default for this is every 100 steps. The output group generic/steps contains an array to which actual step of the MD trajectory a snapshot belongs.

temperatures = job_loaded['output/generic/temperature']
steps = job_loaded['output/generic/steps']
plt.plot(steps, temperatures)
plt.xlabel('MD step')
plt.ylabel('Temperature [K]');

We can also make a histogram of the positions during the simulation.

pos = job_loaded['output/generic/positions']
# The first axis refers to the simulation step, the second to the atom in the structure.
x = pos[:, :, 0]
y = pos[:, :, 1]
z = pos[:, :, 2]
# Let's pick only the bottom layer of atoms
sel = np.abs(z) < 0.1
plt.hist2d(x[sel], y[sel], bins=50)
plt.xlabel('x [$\AA$]')
plt.ylabel('y [$\AA$]')
plt.gca().set_aspect('equal', 'box');

Running an atomistic calculation using DFT (with SPHInX)

job_sphinx = pr.create.job.Sphinx("sphinx_job")
job_sphinx.structure = pr.create.structure.bulk('Al', a=4.04)
job_sphinx.set_exchange_correlation_functional("PBE")
job_sphinx.plane_wave_cutoff = 350
job_sphinx.set_occupancy_smearing("fermi", width=0.2)
job_sphinx.calc_static()
job_sphinx.run()

job_sphinx["output/generic"]

job_sphinx["output/generic/energy_tot"] # Energy for every ionic step

job_sphinx['output/electronic_structure']

job_sphinx.input

eigenvalues = job_sphinx['output/electronic_structure/eig_matrix'].flatten()
occupancies = job_sphinx['output/electronic_structure/occ_matrix'].flatten()

efermi = job_sphinx['output/electronic_structure/efermi']
args = np.argsort(eigenvalues)
# args = np.argsort(eigenvalues)
plt.plot(eigenvalues[args] - efermi, occupancies[args], "-x")
plt.xlabel('Energy')
plt.ylabel('Occupancy');

Additional Reading

Previous workshop materials:

• Phasediagram Workshop 2020

• Potentials Workshop 2021

• CECAM Workshop 2021

• Potentials Workshop 2022

• Youtube

Software from us:

• pyiron github

• eisenforschung github

• pyiron for experimental data

• pyiron for continuum simulations

• Automatic Melting Point Calculations

• pyiron for machine learning potentials

For Take-Away

More Structure Creation Tools

Creating a bulk fcc primitive unitcell and supercell.

Al_unitcell_primitive = pr.create.structure.bulk('Al', a=4.04)
Al_unitcell_primitive.repeat([3, 3, 3]).plot3d(particle_size=2)

Creating a vacancy is easy, just delete an atom.

Al_vacancy = Al_supercell_3_3_3.copy()
del Al_vacancy[0] # Deleting the first atom
print(Al_supercell_3_3_3.get_chemical_formula(), Al_vacancy.get_chemical_formula())
Al_vacancy.plot3d(particle_size=2)

Surfaces are also easy to create, thanks to ASE integration.

num_layers = 4
Al_fcc_111 = pr.create.structure.surface("Al", surface_type="fcc111", size=(4,4,num_layers), vacuum=12, orthogonal=True)
Al_fcc_111.plot3d(particle_size=2)

We also have some specialized methods for compounds.

pr.create.structure.compound.B2("Li", "Al").repeat(2).plot3d()

If needed you can also have full control over the structure. Here we’ll use it to setup a single atom and a dimer.

cell = np.eye(3) * 10
Al_atom_box = pr.create.structure.atoms("Al", cell=cell, scaled_positions=[[0.5, 0.5, 0.5]])
Al_atom_box.plot3d(particle_size=2)

cell = np.eye(3) * 10
Al_atom_1 = pr.create.structure.atoms("Al", cell=cell, scaled_positions=[[0.5, 0.5, 0.5]])
Al_atom_2 = Al_atom_1.copy()

dimer_length = 2.5
Al_atom_2.positions[:, 2] += dimer_length
Al_dimer = Al_atom_1 + Al_atom_2
Al_dimer.center()
Al_dimer.plot3d()

This demonstrates also how to build more complex structures, just + them together, if they have the same boundary conditions.

Energy volume curve for Al

First we create some functions to extract the final energy and volume of a job

def get_volume(job):
    return job["output/generic/volume"][-1]

def get_energy(job):
    return job["output/generic/energy_tot"][-1]

job_loaded = pr['lammps_job']

get_energy(job_loaded), get_volume(job_loaded)

We create first a few jobs with different lattice parameters.

Project can be nested within each other with create_group.

pr_ev = pr.create_group("E_V_curve") # Creating a new sub-project within the main project
a_list = np.linspace(3.8, 4.4, 7)
for a in a_list:
    job_name = "job_a_{:.4}".format(a).replace(".", "_")
    job = pr_ev.create.job.Lammps(job_name)
    job.structure = pr_ev.create.structure.bulk("Al", a=a)
    job.potential = "2005--Mendelev-M-I--Al-Fe--LAMMPS--ipr1"
    job.calc_minimize()
    job.run()

pr

pr_ev.job_table()

Now that we have all the jobs and they finished correctly, we can analyze them.

vol_list = []
energy_list = []

for job in pr["E_V_curve"].iter_jobs(status="finished"):
    vol_list.append(get_volume(job))
    energy_list.append(get_energy(job))

args = np.argsort(vol_list)
vol_list = np.array(vol_list)
energy_list = np.array(energy_list)
plt.plot(vol_list[args], energy_list[args], "-x")
plt.xlabel("Volume [$\mathrm{\AA^3}$]")
plt.ylabel("Energy [eV]");

E-V curves for DFT

We now do the same with a DFT job to allow us to compare. Notice that the structure of the code is almost the same between DFT and Lammps, just the input parameters are different.

pr_ev = pr.create_group("E_V_curve_DFT") # Creating a new sub-project within the main project
a_list = np.linspace(3.8, 4.4, 7)
for a in a_list:
    job_name = "job_a_{:.4}".format(a).replace(".", "_")
    job = pr_ev.create.job.Sphinx(job_name)
    job.structure = pr_ev.create.structure.bulk("Al", a=a)
    job.set_exchange_correlation_functional("PBE")
    job.plane_wave_cutoff = 350
    job.set_kpoints(mesh=[4, 4, 4])
    job.set_occupancy_smearing("fermi", 0.1)
    job.calc_minimize()
    job.run()

Also for the analysis the code necessary is the same, allowing us to easily loop over both groups and plot them together.

for group in pr.iter_groups():
    vol_list = list()
    energy_list = list()
    for job in group.iter_jobs(status="finished"):
        vol_list.append(get_volume(job))
        energy_list.append(get_energy(job))
    args = np.argsort(vol_list)
    vol_list = np.array(vol_list)
    energy_list = np.array(energy_list)
    plt.plot(vol_list[args], energy_list[args] - np.min(energy_list), "-x", label=group.name)
plt.xlabel("Volume [$\mathrm{\AA^3}$]")
plt.ylabel("Energy [eV]")
plt.legend();

Advanced pyiron: Automated workflows and analysis tools

While we could in principle obtain thee E-V cureves by setting up and analyzing the calculations manually as done above, we could also use predefined workflows in pyiron which does this automatically

num_pot = 3
pot_finder = pr.inspect_emperical_potentials()
pot_list = pot_finder.find("Al").Name.to_list()[20:20+num_pot]
pot_list

def clean_project_name(name):
    return name.replace('-', '_')

# Automated Murnaghan using pyiron tables

pr_murn = Project("murn_auto") 

for pot in pot_list:
    
    lammps_job = pr_murn.create.job.Lammps("lammps_ref")
    lammps_job.structure = pr.create.structure.bulk("Al")
    lammps_job.potential = pot
    lammps_job.calc_minimize()
    
    # Creating a Murnaghan workflow (char names not to exceed 50 chars)
    job_name = f"murn_ref_{clean_project_name(pot)}"[:40]
    
    # The job type 'Murnaghan' sets up the appropriate workflow 
    murn_job = lammps_job.create_job(pr.job_type.Murnaghan, job_name)
    murn_job.input["num_points"] = 9
    murn_job.run()

murn_job.plot()

murn_job["output/equilibrium_volume"], murn_job["output/equilibrium_bulk_modulus"]

np.linalg.norm(murn_job["output/structure/cell/cell"][0]) * np.sqrt(2)

We now analyze the data using our in-built pyiron tables class

# A filter function that selects only Murnaghan jobs
def get_only_murn(job_table):
    return (job_table.hamilton == "Murnaghan") & (job_table.status == "finished")

# Functions to obtain output from Murnaghan jobs
def get_eq_vol(job_path):
    return job_path["output/equilibrium_volume"]

def get_eq_lp(job_path):
    return np.linalg.norm(job_path["output/structure/cell/cell"][0]) * np.sqrt(2)

def get_eq_bm(job_path):
    return job_path["output/equilibrium_bulk_modulus"]

def get_potential(job_path):
    return job_path["lammps_ref/input/potential/Name"]

# Creating a pyiron table and processing output
table = pr_murn.create.table("table_murn", delete_existing_job=True)
table.db_filter_function = get_only_murn
table.add["potential"] = get_potential
table.add["a"] = get_eq_lp
table.add["eq_vol"] = get_eq_vol
table.add["eq_bm"] = get_eq_bm

table.run()
table.get_dataframe()

Extra Credits

1. Check the E-V curves for different crystal structures of the same element with Murnaghan jobs (Hint: execute pr.create.structure.bulk? in a cell)

2. Calculate the thermal expansion for an element and structure of your liking (Hint: run MD at different temperatures with the argument pressure=0)