Exercise 1: 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.

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 ‘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

'/home/jovyan/potentials/introduction/first_steps/'

pr

{'groups': [], 'nodes': []}

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

Al: [0. 0. 0.]
Al: [0.   2.02 2.02]
Al: [2.02 0.   2.02]
Al: [2.02 2.02 0.  ]
pbc: [ True  True  True]
cell: 
Cell([4.04, 4.04, 4.04])

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 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)

Al108 Al107

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.

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]

['2004--Mishin-Y--Ni-Al--LAMMPS--ipr2',
 '2004--Zhou-X-W--Al--LAMMPS--ipr2',
 '2005--Mendelev-M-I--Al-Fe--LAMMPS--ipr1',
 '2007--Silva-A-C--Al-Ni--LAMMPS--ipr1',
 '2008--Mendelev-M-I--Al--LAMMPS--ipr1',
 '2009--Kim-Y-M--Mg-Al--LAMMPS--ipr1',
 '2009--Mendelev-M-I--Al-Mg--LAMMPS--ipr1',
 '2009--Purja-Pun-G-P--Ni-Al--LAMMPS--ipr1',
 '2009--Zhakhovskii-V-V--Al--LAMMPS--ipr2',
 '2010--Lee-E--Fe-Al--LAMMPS--ipr1']

# 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
Parameter Value Comment
0 units metal
1 dimension 3
2 boundary p p p
3 atom_style atomic
4 read_data structure.inp
5 include potential.inp
6 fix___ensemble all npt temp 800.0 800.0 0.1 iso 0.0 0.0 1.0
7 variable___dumptime equal 100
8 variable___thermotime equal 100
9 timestep 0.001
10 velocity all create 1600.0 13455 dist gaussian
11 dump___1 all custom ${dumptime} dump.out id type xsu ysu zsu fx fy fz vx vy vz
12 dump_modify___1 sort id format line "%d %d %20.15g %20.15g %20.15g %20.15g %20.15g %20.15g %20.15g %20.15g %20.15g"
13 thermo_style custom step temp pe etotal pxx pxy pxz pyy pyz pzz vol
14 thermo_modify format float %20.15g
15 thermo ${thermotime}
16 run 10000

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()

The job lammps_job was saved and received the ID: 47

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

pr

{'groups': [], 'nodes': ['lammps_job']}

You can get a quick overview with the job_table method.

pr.job_table()
id status chemicalformula job subjob projectpath project timestart timestop totalcputime computer hamilton hamversion parentid masterid
0 47 finished Al108 lammps_job /lammps_job /home/jovyan/ potentials/introduction/first_steps/ 2022-06-07 16:22:34.601906 2022-06-07 16:22:38.109900 3.0 pyiron@jupyter-m-2epoul#1 Lammps 0.1 None None

Analysing a calculation#

# Load the job
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

{'groups': ['executable', 'input', 'output'], 'nodes': ['HDF_VERSION', 'NAME', 'TYPE', 'VERSION', 'server', 'status']}

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

job_loaded["output"]

{'groups': ['generic', 'structure'], 'nodes': []}

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"]

{'groups': [], 'nodes': ['cells', 'energy_pot', 'energy_tot', 'forces', 'indices', 'positions', 'pressures', 'steps', 'temperature', 'unwrapped_positions', 'velocities', 'volume']}

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

[<matplotlib.lines.Line2D at 0x7f101a8870a0>]
../_images/ex_01_introduction_to_pyiron_56_1.png

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]');
../_images/ex_01_introduction_to_pyiron_62_0.png

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');
../_images/ex_01_introduction_to_pyiron_64_0.png

Running an atomistic calculation using DFT (with SPHInX)#

job_sphinx = pr.create.job.Sphinx("sphinx_job")
job_sphinx.structure = Al_unitcell_primitive
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()

The job sphinx_job was saved and received the ID: 48

job_sphinx["output/generic"]

DataContainer({
  "dft": {
    "n_valence": "{'Al': 3}",
    "bands_k_weights": "array([0.03125, 0.09375, 0.09375, 0.09375, 0.09375, 0.1875 , 0.1875 ,\n       0.09375, 0.03125, 0.09375])",
    "kpoints_cartesian": "array([[ 0.15999579,  0.15999579,  0.15999579],\n       [-0.15999579, -0.15999579,  1.11997052],\n       [-0.47998736, -0.47998736,  2.07994525],\n       [-0.79997894, -0.79997894,  3.03991998],\n       [-0.47998736,  0.79997894,  0.79997894],\n       [-0.79997894,  0.47998736,  1.75995367],\n       [-1.11997052,  0.15999579,  2.7199284 ],\n       [-1.11997052,  1.43996209,  1.43996209],\n       [ 0.47998736,  0.47998736,  0.47998736],\n       [ 0.15999579,  0.15999579,  1.43996209]])",
    "bands_e_fermi": "array([5.826281, 5.728571])",
    "bands_occ": "array([[[2.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00,\n         0.0000e+00],\n        [2.0000e+00, 2.0000e-04, 0.0000e+00, 0.0000e+00, 0.0000e+00,\n         0.0000e+00],\n        [2.0000e+00, 2.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00,\n         0.0000e+00],\n        [2.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00,\n         0.0000e+00],\n        [2.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00,\n         0.0000e+00],\n        [2.0000e+00, 2.0172e+00, 5.3000e-03, 0.0000e+00, 0.0000e+00,\n         0.0000e+00],\n        [2.0000e+00, 8.9180e-01, 0.0000e+00, 0.0000e+00, 0.0000e+00,\n         0.0000e+00],\n        [2.0000e+00, 1.8580e-01, 2.0000e-04, 0.0000e+00, 0.0000e+00,\n         0.0000e+00],\n        [2.0000e+00, 2.0212e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00,\n         0.0000e+00],\n        [2.0000e+00, 2.0002e+00, 1.0000e-04, 0.0000e+00, 0.0000e+00,\n         0.0000e+00]]])",
    "bands_eigen_values": "array([[[ 5.1303, 14.0305, 17.448 , 17.448 , 21.4033, 21.4038],\n        [ 2.886 ,  7.8761, 12.2759, 17.483 , 19.4736, 19.566 ],\n        [ 0.4173,  2.6042, 13.2754, 14.5396, 16.3723, 20.5521],\n        [ 4.0041, 10.9748, 14.5241, 16.1197, 17.3349, 19.8625],\n        [ 1.7777,  8.874 , 11.2518, 12.5463, 14.1234, 18.1266],\n        [ 2.5362,  4.6844,  7.0896,  9.2511, 17.583 , 21.7615],\n        [ 0.6753,  5.7575, 10.227 , 12.3039, 18.524 , 20.6616],\n        [ 1.4632,  6.0015,  7.8495, 11.5583, 13.0563, 17.1834],\n        [ 1.7756,  4.7487, 16.5213, 16.5215, 18.5486, 18.5497],\n        [ 1.5011,  3.6299,  8.1685, 13.5452, 15.4597, 23.5714]]])",
    "scf_convergence": "[True]",
    "scf_energy_int": "[[-57.07249516929677, -57.186612261850335, -57.187335449229764, -57.187341601887454]]",
    "scf_energy_free": "[[-57.0952651989044, -57.209309192558685, -57.21005766552014, -57.21006554754307]]",
    "scf_magnetic_forces": "[]",
    "scf_computation_time": "[[0.264505, 0.528669, 0.827663, 1.125228]]",
    "scf_energy_zero": "[[-57.08957269148888, -57.20363495987479, -57.20437711144075, -57.20438456112916]]",
    "scf_energy_band": "[[3.74461494688908, 3.5621788078508088, 3.576114113284531, 3.5755812775248286]]",
    "scf_electronic_entropy": "[[3.098019981452102, 3.088074371020078, 3.0915146481058895, 3.091749941969506]]",
    "scf_residue": "[[0.00499972, 0.000815417, 5.41578e-05, 1.06328e-05]]",
    "energy_int": "array([-57.1873416])",
    "energy_free": "array([-57.21006555])",
    "magnetic_forces": "array([], dtype=float64)",
    "computation_time": "array([1.125228])",
    "energy_zero": "array([-57.20438456])",
    "energy_band": "array([3.57558128])",
    "electronic_entropy": "array([3.09174994])",
    "residue": "array([1.06328e-05])"
  },
  "positions": "array([[[0., 0., 0.]]])",
  "cells": "array([[[0.  , 2.02, 2.02],\n        [2.02, 0.  , 2.02],\n        [2.02, 2.02, 0.  ]]])",
  "volume": "array([16.48480823])",
  "forces": "array([[[-0.,  0.,  0.]]])",
  "energy_tot": "array([-57.21006555])",
  "energy_pot": "array([-57.21006555])"
})
job_sphinx["output/generic/energy_tot"] # Energy for every ionic step

array([-57.21006555])

job_sphinx['output/electronic_structure']

ElectronicStructure Instance
----------------------------
Number of spin channels: 1
Number of k-points: 10
Number of bands: 6
spin 0: Is a metal: True

job_sphinx.input

Group({
  "sphinx": {
    "pawPot": {
      "species": [
        {
          "name": "'\"Al\"'",
          "potType": "'\"AtomPAW\"'",
          "element": "'\"Al\"'",
          "potential": "'\"Al_GGA.atomicdata\"'"
        }
      ]
    },
    "structure": {
      "cell": "array([[0.        , 3.81724677, 3.81724677],\n       [3.81724677, 0.        , 3.81724677],\n       [3.81724677, 3.81724677, 0.        ]])",
      "species": [
        {
          "element": "'\"Al\"'",
          "atom": [
            {
              "coords": "array([0., 0., 0.])"
            }
          ]
        }
      ]
    },
    "basis": {
      "eCut": "25.724525522958494",
      "kPoint": {
        "coords": "array([0.5, 0.5, 0.5])",
        "weight": "1",
        "relative": "True"
      },
      "folding": "array([4, 4, 4])",
      "saveMemory": "True"
    },
    "PAWHamiltonian": {
      "nEmptyStates": "4",
      "ekt": "0.2",
      "FermiDirac": "1",
      "xc": "'PBE'",
      "spinPolarized": "False"
    },
    "initialGuess": {
      "waves": {
        "lcao": [],
        "pawBasis": "True"
      },
      "rho": {
        "atomicOrbitals": "True"
      },
      "noWavesStorage": "False"
    },
    "main": {
      "scfDiag": [
        {
          "rhoMixing": "'1.0'",
          "spinMixing": "'1.0'",
          "dEnergy": "3.674932217565499e-06",
          "maxSteps": "'100'",
          "preconditioner": {
            "type": "'KERKER'",
            "scaling": "1.0",
            "spinScaling": "1.0"
          },
          "blockCCG": []
        }
      ],
      "evalForces": {
        "file": "'\"relaxHist.sx\"'"
      }
    }
  },
  "EnCut": "350",
  "KpointCoords": "[0.5, 0.5, 0.5]",
  "KpointFolding": "[4, 4, 4]",
  "EmptyStates": "4",
  "Sigma": "0.2",
  "Xcorr": "'PBE'",
  "VaspPot": "False",
  "Estep": "100",
  "Ediff": "0.0001",
  "WriteWaves": "True",
  "KJxc": "False",
  "SaveMemory": "True",
  "rhoMixing": "1.0",
  "spinMixing": "1.0",
  "rhoResidualScaling": "1.0",
  "spinResidualScaling": "1.0",
  "CheckOverlap": "True",
  "THREADS": "1",
  "use_on_the_fly_cg_optimization": "True",
  "FermiDirac": "1"
})
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');
../_images/ex_01_introduction_to_pyiron_71_0.png

Task 1: 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)

(-339.86911482815, 1897.35642816221)

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()

The job job_a_3_8 was saved and received the ID: 49
The job job_a_3_9 was saved and received the ID: 50
The job job_a_4_0 was saved and received the ID: 51
The job job_a_4_1 was saved and received the ID: 52
The job job_a_4_2 was saved and received the ID: 53
The job job_a_4_3 was saved and received the ID: 54
The job job_a_4_4 was saved and received the ID: 55

pr

{'groups': ['E_V_curve'], 'nodes': ['lammps_job', 'sphinx_job']}

pr_ev.job_table()
id status chemicalformula job subjob projectpath project timestart timestop totalcputime computer hamilton hamversion parentid masterid
0 49 finished Al job_a_3_8 /job_a_3_8 /home/jovyan/ potentials/introduction/first_steps/E_V_curve/ 2022-06-07 16:22:47.261636 2022-06-07 16:22:47.860877 0.0 pyiron@jupyter-m-2epoul#1 Lammps 0.1 None None
1 50 finished Al job_a_3_9 /job_a_3_9 /home/jovyan/ potentials/introduction/first_steps/E_V_curve/ 2022-06-07 16:22:49.073172 2022-06-07 16:22:49.352637 0.0 pyiron@jupyter-m-2epoul#1 Lammps 0.1 None None
2 51 finished Al job_a_4_0 /job_a_4_0 /home/jovyan/ potentials/introduction/first_steps/E_V_curve/ 2022-06-07 16:22:50.395300 2022-06-07 16:22:51.172536 0.0 pyiron@jupyter-m-2epoul#1 Lammps 0.1 None None
3 52 finished Al job_a_4_1 /job_a_4_1 /home/jovyan/ potentials/introduction/first_steps/E_V_curve/ 2022-06-07 16:22:51.994665 2022-06-07 16:22:52.415024 0.0 pyiron@jupyter-m-2epoul#1 Lammps 0.1 None None
4 53 finished Al job_a_4_2 /job_a_4_2 /home/jovyan/ potentials/introduction/first_steps/E_V_curve/ 2022-06-07 16:22:53.700368 2022-06-07 16:22:54.143684 0.0 pyiron@jupyter-m-2epoul#1 Lammps 0.1 None None
5 54 finished Al job_a_4_3 /job_a_4_3 /home/jovyan/ potentials/introduction/first_steps/E_V_curve/ 2022-06-07 16:22:55.563556 2022-06-07 16:22:56.456825 0.0 pyiron@jupyter-m-2epoul#1 Lammps 0.1 None None
6 55 finished Al job_a_4_4 /job_a_4_4 /home/jovyan/ potentials/introduction/first_steps/E_V_curve/ 2022-06-07 16:22:57.806968 2022-06-07 16:22:58.085522 0.0 pyiron@jupyter-m-2epoul#1 Lammps 0.1 None None

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]");
../_images/ex_01_introduction_to_pyiron_82_1.png

Task 2: 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()

The job job_a_3_8 was saved and received the ID: 56
The job job_a_3_9 was saved and received the ID: 57
The job job_a_4_0 was saved and received the ID: 58
The job job_a_4_1 was saved and received the ID: 59
The job job_a_4_2 was saved and received the ID: 60
The job job_a_4_3 was saved and received the ID: 61
The job job_a_4_4 was saved and received the ID: 62

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();
../_images/ex_01_introduction_to_pyiron_86_3.png

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

['2004--Mishin-Y--Ni-Al--LAMMPS--ipr2',
 '2004--Zhou-X-W--Al--LAMMPS--ipr2',
 '2005--Mendelev-M-I--Al-Fe--LAMMPS--ipr1']

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()

The job murn_ref_2004__Mishin_Y__Ni_Al__LAMMPS__ was saved and received the ID: 63
The job murn_ref_2004__Mishin_Y__Ni_Al__LAMMPS___0_9 was saved and received the ID: 64
The job murn_ref_2004__Mishin_Y__Ni_Al__LAMMPS___0_925 was saved and received the ID: 65
The job murn_ref_2004__Mishin_Y__Ni_Al__LAMMPS___0_95 was saved and received the ID: 66
The job murn_ref_2004__Mishin_Y__Ni_Al__LAMMPS___0_975 was saved and received the ID: 67
The job murn_ref_2004__Mishin_Y__Ni_Al__LAMMPS___1_0 was saved and received the ID: 68
The job murn_ref_2004__Mishin_Y__Ni_Al__LAMMPS___1_025 was saved and received the ID: 69
The job murn_ref_2004__Mishin_Y__Ni_Al__LAMMPS___1_05 was saved and received the ID: 70
The job murn_ref_2004__Mishin_Y__Ni_Al__LAMMPS___1_075 was saved and received the ID: 71
The job murn_ref_2004__Mishin_Y__Ni_Al__LAMMPS___1_1 was saved and received the ID: 72
The job murn_ref_2004__Zhou_X_W__Al__LAMMPS__ipr was saved and received the ID: 73
The job murn_ref_2004__Zhou_X_W__Al__LAMMPS__ipr_0_9 was saved and received the ID: 74
The job murn_ref_2004__Zhou_X_W__Al__LAMMPS__ipr_0_925 was saved and received the ID: 75
The job murn_ref_2004__Zhou_X_W__Al__LAMMPS__ipr_0_95 was saved and received the ID: 76
The job murn_ref_2004__Zhou_X_W__Al__LAMMPS__ipr_0_975 was saved and received the ID: 77
The job murn_ref_2004__Zhou_X_W__Al__LAMMPS__ipr_1_0 was saved and received the ID: 78
The job murn_ref_2004__Zhou_X_W__Al__LAMMPS__ipr_1_025 was saved and received the ID: 79
The job murn_ref_2004__Zhou_X_W__Al__LAMMPS__ipr_1_05 was saved and received the ID: 80
The job murn_ref_2004__Zhou_X_W__Al__LAMMPS__ipr_1_075 was saved and received the ID: 81
The job murn_ref_2004__Zhou_X_W__Al__LAMMPS__ipr_1_1 was saved and received the ID: 82
The job murn_ref_2005__Mendelev_M_I__Al_Fe__LAMM was saved and received the ID: 83
The job murn_ref_2005__Mendelev_M_I__Al_Fe__LAMM_0_9 was saved and received the ID: 84
The job murn_ref_2005__Mendelev_M_I__Al_Fe__LAMM_0_925 was saved and received the ID: 85
The job murn_ref_2005__Mendelev_M_I__Al_Fe__LAMM_0_95 was saved and received the ID: 86
The job murn_ref_2005__Mendelev_M_I__Al_Fe__LAMM_0_975 was saved and received the ID: 87
The job murn_ref_2005__Mendelev_M_I__Al_Fe__LAMM_1_0 was saved and received the ID: 88
The job murn_ref_2005__Mendelev_M_I__Al_Fe__LAMM_1_025 was saved and received the ID: 89
The job murn_ref_2005__Mendelev_M_I__Al_Fe__LAMM_1_05 was saved and received the ID: 90
The job murn_ref_2005__Mendelev_M_I__Al_Fe__LAMM_1_075 was saved and received the ID: 91
The job murn_ref_2005__Mendelev_M_I__Al_Fe__LAMM_1_1 was saved and received the ID: 92

murn_job.plot()
../_images/ex_01_introduction_to_pyiron_91_0.png 
<AxesSubplot:title={'center':'Murnaghan: error: 1.1881350053920211e-07'}, xlabel='Volume ($\\AA^3$)', ylabel='energy (eV)'>

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

(16.41961401289247, 85.21427045715596)

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

4.050000000000001

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()

The job table_murn was saved and received the ID: 93
job_id potential a eq_vol eq_bm
0 63 2004--Mishin-Y--Ni-Al--LAMMPS--ipr2 4.05 16.606859 80.886124
1 73 2004--Zhou-X-W--Al--LAMMPS--ipr2 4.05 16.611633 71.503293
2 83 2005--Mendelev-M-I--Al-Fe--LAMMPS--ipr1 4.05 16.419614 85.214270

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)