# 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.

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
job_md.input.control["run"] = ____
```

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
```
```# 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
```
```# 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.