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
Importing necessary libraries #
As a first step we import the libraries numpy for data analysis and matplotlib for visualization.
import numpy as np
%matplotlib inline
import matplotlib.pylab as plt
Fundamentally, we only need to import one module from pyiron: the Project class
from pyiron import Project
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 ‘tutorial’.
Working with atomistic structures #
Creation of a project instance #
pr = Project("tutorial")
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(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)
Creating random crystals #
pyxtal is a program to generate random structures of a given space group and stoichiometry.
It is useful for crystal structure prediction and also for generating training structures.
We will use it later again, so let’s briefly look at how it works.
import structuretoolkit as stk
from pyiron_atomistics import ase_to_pyiron
from pyiron_atomistics.atomistics.structure.structurestorage import StructureStorage
structuretoolkit is a library for structure manipulation and analysis also developed by the pyiron team. For compatibility with a wider range of codes it operates purely on ASE Atoms objects, so we need to convert structures explicitely here. In the next release of pyiron_atomistics you will be able to call pr.create.structure.pyxtal directly for a more convenient wrapper.
fcc = ase_to_pyiron(stk.build.pyxtal(
group=225, # fcc
species=['Al'], # list of atoms we want to place
num_ions=[4], # how of each of them
))
fcc
fcc.plot3d()
groups = [1, 194, 225]
store = StructureStorage()
for structure in stk.build.pyxtal(
group=groups,
species=['Al'],
num_ions=[4]
):
store.add_structure(structure['atoms'])
structure = store.get_structure(1)
structure.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("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 = pr.create.structure.bulk('Al', cubic=True).repeat(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()
Once it is finished we can access the parsed output.
job_lammps['output']
plt.plot(job_lammps['output/generic/energy_pot'])
Collecting structures and energies for training #
To train a machine learning potential one usually starts with a set of structures, then run static DFT on each of them. We will simulate this here, by running lammps with an existing potential to save some time.
Let’s start by adding structures to our container.
store = StructureStorage()
A usual approach is to start with some prototypical structures and strain and shake them. We can use the fcc Al we created earlier.
fcc
strains = np.linspace(-0.05, 0.05, 5)
strains
for eps in strains:
store.add_structure(
structure=fcc.apply_strain(eps, return_box=True),
identifier=f'fcc_{eps}'
)
fcc.rattle?
for eps in strains:
structure = fcc.apply_strain(eps, return_box=True)
structure.rattle(0.05)
print('Strain', eps, 'Volume', structure.get_volume(per_atom=True))
store.add_structure(
structure=structure,
identifier=f'fcc_{eps}_rattle'
)
We usually also use random crystals generated by sampling random space groups. One would use more spacegroups and differently sized unit cells, but we’ll keep it simple here. Details can be found in this paper.
https://journals.aps.org/prb/abstract/10.1103/PhysRevB.107.104103
for i, structure in enumerate(stk.build.pyxtal(
group=groups,
species=['Al'],
num_ions=[4],
tm='metallic'
)):
store.add_structure(
structure=structure['atoms'],
identifier=f'random_{i}'
)
Now we can run our “expensive DFT”.
store.number_of_structures
for i, structure in enumerate(store.iter_structures()):
# pyiron has some opinions of what is a proper "job name" for technical reasons, but
# you can always use this function to translate your favorite one into a "proper" one
name = pr.create.job_name(store['identifier', i])
job = pr.create.job.Lammps(name)
job.potential = "2005--Mendelev-M-I--Al-Fe--LAMMPS--ipr1"
job.structure = structure
job.calc_static()
job.run()
Once our reference calculations are finished, we can collect the results in a container for convenience.
train = pr.create.job.TrainingContainer("Aluminium")
train.include_structure?
train.include_job?
for job in pr.iter_jobs(hamilton='Lammps', status='finished'):
for i in range(job.number_of_structures):
train.include_job(job, iteration_step=i)
“Running” the container simply saves it to disk and to our database. It can also precompute nearest neighbor information.
train.run()
Besides keeping everything in one place, the TrainingContainer also defines a number of plotting functions.
df = train.plot.energy_volume()
train.plot.energy_distance()
plt.xlabel(r'Minimum Nearest Neighbor Distance [$\mathrm{\AA}$]')
plt.ylabel('Energy [eV/atom]')
For fitting potentials within pyiron you can use the container as is, but if you have external tools, you can export the data into a table.
train.to_pandas()