Exercise 2: Simulation Codes

Authors: Jan Janssen, Tilmann Hickel, Jörg Neugebauer (Max-Planck-Institut für Eisenforschung)

The scope of this second exercise is to become familar with:

  • Interatomic potential calculation,

  • the pyiron job management and

  • the aggregation of multiple calculations.

Reminder

In the first session we learned how to create a pyiron project object and then use this pyiron project object to create atomistic structure objects.

# Import the Project object
from pyiron import ____
# Create a Project object instance for a project named atomistics
pr = ____("simulations")
# Create a cubic aluminum fcc structure and repeat it 3 times in each direction
al_fcc = pr.create_ase_bulk(____, _____=True)
al_fcc_repeated = al_fcc.repeat(____)
# Confirm the final structure has 108 atoms by calculating the length of the structure object
____(al_fcc_repeated) == 108

LAMMPS Calculation

The Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) code is used inside pyiron for atomistic simualtion with interatomic potentials. These interatomic potentials approximate the interaction of atoms and can be either fitted to density functional theory (DFT) or experimental results. Still in contrast to density functional theory which scales cubically with the number of atoms interatomic potentials scale linearly with the number of atoms. Therefore we are going to use primarly interatomic potentials in this workshop but most calculations could be executed with a DFT codes as well.

Molecular dynamics calculation

We start with a first molecular dynamic calculation at an ensemble with constant number of atoms, constant volume and contant fixed temperature.

# Create a LAMMPS job object with the job name lmp
job_md = pr.create_job(
    job_type=pr.job_type.Lammps, 
    job_name=____
)
# Assign the fcc aluminium structure to the LAMMPS job object
job_md.structure = _________
# Define an ensemble with constant number of atoms, 
# constant volume and a constant temperature of 500K
# and simulate 10000 molecular dynamics steps
job_md.calc_md(temperature=______, n_ionic_steps=______)
# Execute the calculation - You get a warning that 
# not potential was set and the default potential was
# used instead. The selection of interatomic potentials 
# is discussed below.
job_md.run()
# Animate the molecular dynamics trajectory
job_md.animate_structure()

Reminder: Job Management in pyiron

After the successful execution of the calculation it is listed in the project database and can be reloaded using either the job name or the job id, therefore theese have to be unique for a given project.

# list all calculations in the current project using the project job table
______.job_table()

As discussed in the previous section by default pyiron calculations are reloaded from the database when a calculation with the same name already exists in a given project. Therefore to overwrite the calculation parameters we use delete_existing_job parameter in the run() function.

# Change the temperature to 800K and calculate 20000 steps
# You get a warning message that the job is already finished. 
job_md.calc_md(temperature=_____, n_ionic_steps=_______)
# Execute the LAMMPS calculation by calling the run function 
# with the delete_existing parameter set to true: 
job_md.run(delete_existing_job=______)

Reminder: Plot calculation results using matplotlib

In the same way we can again use matplotlib to analyse the results of our calculation. For example we can use the matplotlib library to plot the temperature over simulation steps.

# Import the matplotlib library for plotting. 
import matplotlib.pyplot as plt

# for the LAMMPS job object plot the temperature over simualation steps
plt.plot(________.output.steps, _____.output.temperature)
plt.xlabel("timesteps")
plt.ylabel("Temperature");

In the beginning the potential energy is close to the 0K equilibrium, therefore to accelerate the equilibration the kinetic energy is set to twice the expected kinetic energy. The additional kinetic energy is transfered to the potential energy resulting in an equal distribution of potential and kinetic energy. With this trick the equilibration is accelerated. The large temperature fluctuations are related to the small number of atoms in the simulation cell.

Advanced input options

Besides the general functions calc_static(), calc_minimize() and calc_md(). pyiron also has the option to modify the input of the simulation code directly.

# display the LAMMPS input file 
job_md.input.control
# change the number of simulation steps to 2000 
# by manually modifying the run command.
# You again get a warning message that the job 
# is already finished. 
job_md.input.control["run"] = ____

Advanced output options

Besides the output properties it is also possible to access the output of a calculation directly from the data interface which is based om the hierachical file format (HDF5) pyiron is using to store the simulation data.

# print content of the job object 
job_md
# print content of the output of the job object
# use strings to specify the path for the data interface 
job_md[_______]
# print content of the generic group 
# of the output of the job object 
job_md[______/______]
# plot the temperature over simulation steps directly from the HDF5 file 
plt.plot(job_md[___________], job_md[___________])
plt.xlabel("timesteps")
plt.ylabel("Temperature");

Accessing the original output files of the LAMMPS code

While pyiron parses most of the output of the simulation codes some users might have the need to access addtional output parameters.

# Decompress the LAMMPS job
_______.decompress()
# Read the LAMMPS output file of the LAMMPS job
_______["log.lammps"]

Beyond a single LAMMPS calculation

While for individual LAMMPS calculation an integrated solution like pyiron is not required, pyiron really shines when it comes to combining multiple calculations. So in the following we iterate over a database of existing interatomic potentials and determine the lattice structure by minimizing the simulation cell.

Filter Interatomic Potential Database

We start by identifying suitable interatomic potentials. By default pyiron already filters the interatomic potentials to only list those which include the required elements. Still the user can further filter the list of available potentials for a given project.

# List interatomic potentials from the NIST repository: 
# https://www.ctcms.nist.gov/potentials/
# which include interactions for aluminium by calling view_potentials() on the LAMMPS job object
potential_df = ______.view_potentials()
potential_df[potential_df.Model == "NISTiprpy"]

Choose an interatmoic potential

Interatomic potentials can be fitted for specific applications, therefore before selecting a given interatomic potential it is recommended to test basic properties of the potential. In this example we calculate the 0K equilibrium lattice constant by optimizing the atmoic supercell.

Job Template

To develop a scalable simulation protocol, we first define a job template and then apply this template to the available interatomic potentials.

# To compare different potentials we start by creating a template job named lmp_template
job_template = pr.create_job(
    job_type=pr.job_type.Lammps, 
    job_name=______
)
# We assign a cubic fcc aluminium supercell 
job_template.structure = pr.create_ase_bulk(_____, _____)
# Enable volume minimization by specifying the pressure as zero
job_template.calc_minimize(pressure=_____)

Iterate over interatomic potentials

After the template is constructed we can iterate over the database of existing interatomic potentials. In this example we limit the total umber of potentials to three to accelerate the calculations.

# Select the first three potentials from the NIST database
potential_df[potential_df.Model == "NISTiprpy"].Name.values[_____]
# We then iterate over the first three potentials of the NIST database
for p in potential_df[potential_df.Model == "NISTiprpy"].Name.values[_____]:
    # create a copy for each of the template job for each of the potentials
    # without creating a new database entry by setting the new_database_entry 
    # parameter to False 
    job_minimize = job_template.copy_to(
        new_job_name="lmp_" + p.replace("-", "_"), 
        new_database_entry=_____
    )
    
    # We then assign the potential 
    job_minimize.potential = p
    
    # Execute the calculation and set the delete_existing_job parameter to True
    job_minimize.run(delete_existing_job=____)
# list all calculations in the current project using the project job table
pr.job_table()

Validate calculations

We validate the simulation results by confirming the calculations have been executed successfully. In this example one calculation failed.

# We use the pyiron job table of the project object 
# to validate all jobs finished successfully with 
# the status "finished"
df = ____.job_table()
len(df[df.status == _______]) == 4

Analyse results

Iterate over the successful calculation and compare the calculated lattice constants.

# load a job object by the job name in the column job 
# of the pyiron job table 
job = pr.load(______)
# print the final lattice constant 
job["output/generic/cells"][______]
# We iterate over the jobs in a the current project 
for job in pr.iter_jobs():
    # Filter job using only those who have the job status 
    # finished and "lmp_" in the job_name
    if job.status == _____ and ____ in job.job_name: 
        # Print the job name and the lattice constant 
        print(job.job_name, job["output/generic/cells"][_____])

pyiron table

To automate the collection of data from individual calculations pyiron includes the pyiron table object.

# Create a pyiron table object
table = pr.create_table()
# Implement a filter function, which returns true 
# for finished jobs and jobs with "lmp_" in the job_name
def filter_jobs(job):
    return job.status == ______ and ______ in job.job_name
# Implement an analysis functions, which takes a job object
# as an input and returns the lattice constant. Based on the
# previous cell above to print the lattice constant.
def get_lattice_constant(job): 
    return _________________
# Implement a second analysis functions, which takes a 
# job object as an input and returns the job_name.
def get_job_name(job):
    return _________________
# Assigne the filter functions and the analysis functions 
# to the pyiron table object
table.filter_function = filter_jobs
table.add["job_name"] = get_job_name
table.add["lattice_constant"] = get_lattice_constant
# Execute the pyiron table just like a pyiron job object
table.run()
# Return a pandas DataFrame with the collected results 
table.get_dataframe()

Summary

In this section you learned:

  • to execute LAMMPS calculation with generic and code-specifc input,

  • to iterate over multiple interatomic potentials

  • and aggregate the data in a pyiron table object. This technique was afterwards used to calculate the lattice constant of multiple interatomic potentials.

Additional exercises:

  • Calculate the lattice constant at a finite temperature of 500K by averaging the lattice constant over 1000 time steps.

  • Calculate the thermal expansion for multiple potentials.

  • How does an interstitial element change the thermal expansion? How does a vacancy change the thermal expansion?