Main Computational Code
The next step is to write the main computational code. It must be contained in a Python file named test_driver/test_driver.py,
although you may include as many other Python files as you wish for utility functions.
Note
Your test_driver directory should be structured as a Python package. This means that it contains a (probably empty) __init__.py,
and you should use relative imports to import any other Python files from your test_driver.py. For example, if you want to put
some function my_helper_function in a Python file helper_functions.py, you shoud structure your test_driver directory like this:
test_driver/
├─ __init__.py
├─ test_driver.py
├─ helper_functions.py
Then, in your test_driver.py you would import my_helper_function like this:
from .helper_functions import my_helper_function
You must create a class named TestDriver inheriting from SingleCrystalTestDriver.
In your TestDriver class, you must overload the function _calculate().
Besides self, the function must also accept **kwargs. Before **kwargs, you may add any additional arguments that
you wish users or the OpenKIM Pipeline to be able to vary. Try to have default values for as many arguments as possible.
Ideally, only physically meaningful variables without a reasonable default should be mandatory to specify.
Note
Temperature and stress/pressure are commonly used, so they do not need to be added as additional arguments. If needed,
your Test Driver will automatically accept arguments temperature_K and cell_cauchy_stress_eV_angstrom3 or
pressure_eV_angstrom3. Stress and pressure are mutually exclusive -- the pressure argument is simply
a convenient way to specify a hydrostatic stress state. The base class only provides recordkeeping of these
variables, it is up to your code to actually impose these conditions on the simulation. See the
Example test_driver.py written with ASE for an example of how to access them from inside your _calculate function.
An example test_driver.py from https://github.com/openkim-hackathons/CrystalGenomeASEExample__TD_000000654321_000 is shown below. The _calculate function computes the energy vs. volume
curve for isotropic expansion and compression of a crystal at zero temperature. You can use it as a starting point for your Test Driver.
Additional documentation about functionality not demonstrated in the example can be found at the bottom of this page.
Once test_driver.py is written, you will have to make some small modifications to the auxiliary files provided in the example
Test Driver in order for your Driver to work in the OpenKIM pipeline and the kimvv package: Metadata and KIM Processing Pipeline Integration, kimvv Python package integration. Before this, you will likely wish
to debug your code by invoking it directly in Python. See how to do that here: Testing and Debugging Your Crystal Genome Test Driver
Please read the code and the comments carefully, as they explain many important aspects of writing your own Test Driver. You should understand the usage of the following functions. Click the links below for more information on each:
self._get_atoms(): returns anAtomsobject representing a primitive unit cell of the crystal as a starting point for your calculations.self._update_nominal_parameter_values(): if your Test Driver changes the crystal structure, pass this a primitive unit cell of the crystal to update the nominal crystal structure.Note
If the symmetry of your crystal has changed in the process of your simulation, this function will raise an error. This is intended and it is expected that you do not handle this exception. It is our convention that if the crystal undergoes a symmetry-changing phase transition, the result is invalid and the Test Driver should exit with an error.
self._verify_unchanged_symmetry(): You may also find it useful to check for a symmetry change without changing the stored description of the crystal. This is useful, for example, if your Test Driver is changing temperature or stress around some reference value, and you wish to make sure that the state changes have not induced a phase transition. This function takes anAtomsobject and returnsTrueif the symmetry has not changed, otherwiseFalse.self._add_property_instance_and_common_crystal_genome_keys(): Use this to initialize an Instance of the KIM Property(s) you defined in Creating a Crystal Genome Property. It will automatically be populated with the keys describing the nominal state of the crystal.self._add_key_to_current_property_instance(): Use this to add additional keys to the last Property Instance you created.self._get_temperature()self._get_nominal_crystal_structure_npt(): If you need fine-grained access to the symmetry-reduced AFLOW prototype designation of the crystalself._add_file_to_current_property_instance(): For adding keys with the "file" type to your Property Instance, e.g. restart files or images
Note
The units you choose for writing property key values must be comprehensible by the GNU units utility.
Example test_driver.py written with ASE
"""
Example OpenKIM Crystal Genome Test Driver
==========================================
This is an example demonstrating usage of the kim-tools package. See
https://kim-tools.readthedocs.io for more information.
"""
import shutil
from os import path
from kim_tools import (
SingleCrystalTestDriver,
get_stoich_reduced_list_from_prototype,
minimize_wrapper,
)
class TestDriver(SingleCrystalTestDriver):
def _calculate(
self, max_volume_scale: float = 1e-2, num_steps: int = 10, **kwargs
) -> None:
"""
Computes the energy vs. volume curve for isotropic expansion and compression.
Args:
max_volume_scale:
Maximum fractional change in volume to investigate
num_steps:
Number of steps to take in each direction
"""
# The base class provides ``self._get_atoms()``, which provides a starting
# point for your calculations. This Atoms object is a primitive cell in the
# setting defined in https://doi.org/10.1016/j.commatsci.2017.01.017.
#
# Perform your calculations using this copy. If using ASE, the calculator is
# already attached. Note that if you invoke the Test Driver using an Atoms
# object, this is not the same object. It is a re-generated object produced
# from the symmetry analysis of the input Atoms object. Additionally, any
# changes made to this Atoms object will not be automatically reflected in the
# output of this Test Driver. If your Test Driver changes the nominal crystal
# structure (through equilibration or relaxation), see how to report that in
# the comments preceding the invocation of
# ``self._update_nominal_parameter_values()`` below.
original_atoms = self._get_atoms()
original_cell = original_atoms.get_cell()
num_atoms = len(original_atoms)
# The base class also provides an instance of the OpenKIM
# ``crystal-structure-npt`` property. This is a symmetry-reduced description of
# the nominal crystal structure. See the definition of this property at
# https://openkim.org/properties/show/2023-02-21/staff@noreply.openkim.org/crystal-structure-npt
# For a general explanation of the structure of a KIM Property Instance, see
# https://openkim.org/doc/schema/properties-framework/
# To access a copy of the property instance as a Python dictionary, use
# ``self._get_nominal_crystal_structure_npt()```
# Here we will use the prototype label to infer the number of atoms per
# stoichiometric formula
prototype_label = self._get_nominal_crystal_structure_npt()["prototype-label"][
"source-value"
]
num_atoms_in_formula = sum(
get_stoich_reduced_list_from_prototype(prototype_label)
)
# Because temperature and stress/pressure are such commonly used parameters, we
# provide utility functions to access them. This example Test Driver does not
# use them, this is for demonstration. The units requested here are the
# defaults, they will be what you get if you omit the ``unit`` argument.
# The base class only provides recordkeeping of these
# variables, it is up to your code to actually impose these conditions.
print(f'Temperature (K): {self._get_temperature(unit="K")}')
print(
"Stress (eV/angstrom^3): "
+ str(self._get_cell_cauchy_stress(unit="eV/angstrom^3"))
)
# Internally, only a stress state is tracked, and ``_get_pressure`` gets
# the pressure from the stress state. If `enforce_hydrostatic` is
# `True` (default), an error will be raised if the stress state is not
# hydrostatic.
print(
"Pressure (eV/angstrom^3): "
+ str(self._get_pressure(unit="eV/angstrom^3", enforce_hydrostatic=True))
)
binding_potential_energy_per_atom = []
binding_potential_energy_per_formula = []
volume_per_atom = []
volume_per_formula = []
step_size = max_volume_scale / num_steps
disclaimer = None
print("\nPerforming energy scan...\n")
for i in range(-num_steps, num_steps + 1):
volume_scale = 1 + step_size * i
linear_scale = volume_scale ** (1 / 3)
# Get atoms again
atoms = self._get_atoms()
atoms.set_cell(original_cell * linear_scale, scale_atoms=True)
volume = atoms.get_volume()
current_volume_per_atom = volume / num_atoms
problem_occurred = False
try:
# The Atoms object returned by ``self._get_atoms()`` comes
# pre-initialized with an ASE calculator. If you need to access the
# name of the KIM model (for example, if you are exporting the atomic
# configuration to run an external simulator like LAMMPS), it can
# be accessed at ``self.kim_model_name``
minimize_wrapper(atoms, variable_cell=False)
# ``self._verify_unchanged_symmetry()`` is used to check that the
# material has not undergone a phase transition.
# It is OpenKIM convention that a run of SingleCrystalTestDriver should
# not involve phase transitions, so here any points where the crystal
# symmetry has changed are not recorded in the curve.
if self._verify_unchanged_symmetry(atoms):
potential_energy = atoms.get_potential_energy()
stress = atoms.get_stress()
print(f"Volume: {volume} Energy: {potential_energy}")
# If your Test Driver changes the nominal crystal structure (e.g.
# through relaxation or MD) AND you intend to write the changed
# structure as a property instance, you must update the nominal
# crystal structure before writing properties. This is done by
# providing an Atoms object to
# ``self._update_nominal_crystal_structure()``.
# SingleCrystalTestDriver expects the crystal to maintain the same
# space group and occupied Wyckoff positions, meaning only the free
# parameters of the crystal unconstrained by symmetry are allowed
# to change. You must provide an Atoms object that is a primitive
# cell in the same setting. Translations, permutations, and rigid
# rotations of the unit cell are permissible, but the identity of
# the lattice vectors may not change during the property computation
# (e.g. lattice vectors a and c may not exchange places, even if
# allowed by symmetry). ``_update_nominal_parameter_values``
# raises an exception if a symmetry change is detected. Here, it
# will be caught and the Test Driver will move on to the next
# structure evaluation. If your Test Driver only deals with one
# structure and it changes symmetry, you should allow it to fail.
self._update_nominal_parameter_values(atoms)
# Normally, you would use either ``_verify_unchanged_symmetry``
# (if you don't intend to write the changed structure), or
# ``_update_nominal_parameter_values`` (if you do), not both, as
# done here for demonstration.
else:
print("Atoms underwent symmetry change")
problem_occurred = True
except Exception:
print(
"Failed to get energy, stress, or parameter values "
f"at volume {volume}"
)
problem_occurred = True
if problem_occurred:
disclaimer = (
"At least one of the requested deformations of the unit cell "
"underwent a phase transformation or otherwise failed to evaluate."
)
else:
current_binding_potential_energy_per_atom = potential_energy / num_atoms
volume_per_atom.append(current_volume_per_atom)
volume_per_formula.append(
current_volume_per_atom * num_atoms_in_formula
)
binding_potential_energy_per_atom.append(
current_binding_potential_energy_per_atom
)
binding_potential_energy_per_formula.append(
current_binding_potential_energy_per_atom * num_atoms_in_formula
)
# Because we have evaluated a structure at some combination of
# temperature and/or stress (only stress in this case), we might as
# well write an instance of the "crystal-structure-npt" property at
# every loop iteration. The newly defined
# "energy-vs-volume-isotropic-crystal" property will be written
# at the end of the loop.
# This method initializes the Property Instance and adds the keys
# common to all Crystal Genome properties. `property_name`` can be the
# full "property-id" field in the Property Definition, or the
# "Property Name", which is just the short name after the slash, as
# used here. The options `write_stress` and `write_temp` can be set to
# `False` (default), `True`, in which case the internally stored nominal
# values will be written, or set to a value in case you wish to write
# a different value. Here we tell the Driver to write a provided stress
# value and write the nominal temperature (zero).
# If you specify the value for stress, you must also specify
# `stress_unit` for the value you have provided.
# For temperature, `temp_unit` defaults to K. See the API documentation
# for the method for more info.
self._add_property_instance_and_common_crystal_genome_keys(
property_name="crystal-structure-npt",
write_stress=stress,
write_temp=True,
stress_unit="eV/angstrom^3",
)
# Because "crystal-structure-npt" has no additional custom fields, we
# are done with this property
# Now it is time to write the output in the format you created in your Property
# Definition. Here we omit stress and temperature because we did not include
# them in our Definition, and optionally add a disclaimer. First, since we have
# been changing the nominal crystal structure, we must reset it to the original
# structure.
self._update_nominal_parameter_values(original_atoms)
self._add_property_instance_and_common_crystal_genome_keys(
property_name="energy-vs-volume-isotropic-crystal",
write_stress=False,
write_temp=False,
disclaimer=disclaimer,
)
# This method adds additional fields to your property instance by specifying the
# key names you defined in your property definition and providing units if
# necessary.
self._add_key_to_current_property_instance(
"volume-per-atom", volume_per_atom, unit="angstrom^3"
)
self._add_key_to_current_property_instance(
"volume-per-formula", volume_per_formula, unit="angstrom^3"
)
# You may also provide a dictionary supplying uncertainty information. It is
# optional, and normally would not be reported for a deterministic calculation
# like this, only one involving some statistics, such as molecular dynamics or
# fitting. There are many possible ways to report uncertainty information,
# detailed at https://openkim.org/doc/schema/properties-framework/
uncertainty_info = {
"source-std-uncert-value": [0] * len(binding_potential_energy_per_atom)
}
self._add_key_to_current_property_instance(
"binding-potential-energy-per-atom",
binding_potential_energy_per_atom,
unit="eV",
uncertainty_info=uncertainty_info,
)
self._add_key_to_current_property_instance(
"binding-potential-energy-per-formula",
binding_potential_energy_per_formula,
unit="eV",
uncertainty_info=uncertainty_info,
)
# This adds a "file" type key to the current Property Instance. It will
# automatically be numbered according to the current Property Instance and
# moved to the output directory if it is not already there
# (e.g. cat.jpg will be automatically moved to the path output/cat-1.jpg
# and reported in the property accordingly.) For this example, we just copy
# the picture packaged with this Test Driver.
cat_path = path.join(path.dirname(path.realpath(__file__)), "data", "cat.jpg")
save_path = "cat.jpg"
shutil.copy(cat_path, save_path)
self._add_file_to_current_property_instance("cat-picture", save_path)
# If your Test Driver reports multiple Property Instances, repeat the process
# above for each one.
Functionality not covered in the above example
self._calc(): This gives you access to the ASE calculator object, if you are building a separateAtomsobject and need to attach a calculator to it.self.kim_model_name: The KIM Model Name (if present). You should only use this if exporting data to a non-ASE simulator (see below).
Writting Test Drivers using LAMMPS
Todo
Write detailed documentation for writing LAMMPS Test Drivers instead of just an overview and example
In general, using ASE to perform the computations is preferable, but you may need access to LAMMPS functionality, for example to run a large parallel MD or static
calculation with domain decomposition. You may use LAMMPS in whichever way is most convenient for you as long as it is wrapped within the _calculate Python
method. For example, one way to run a LAMMPS simulation in this framework is to export your atomic configuration using ase.io.write(), create LAMMPS input file(s)
with kim commands using the KIM model stored in the base class' attribute self.kim_model_name, run your simulation(s), and read the configuration back in using ase.io.read() (for example, to re-detect the changed crystal structure).
If you are running an MD simulation, the structure you report should be time-averaged and averaged over the supercell folded back into the unit cell.
This will give you more robust averages assuming the translational symmetry of the unit cell was not broken. At this time, the function
kim_tools.symmetry_util.core.reduce_and_avg() will perform the averaging assuming your supercell is built from contiguous repeats of the unit cell
(i.e. atoms 0 to N-1 in the supercell are the original unit cell, atoms
N to 2 N-1 are a the original unit cell shifted by an integer multiple of the lattice vectors, and so on). If the translational symmetry
is broken, this function will detect it and raise an error. See https://github.com/openkim/kim-tools/blob/main/tests/test_symmetry_util.py for an example of how to use this function.
Additionally, if you are performing an NPT simulation, you may as well write an instance of the crystal-structure-npt property for future re-use. Note the optional restart-file key. It is recommended that you save a restart file (for example, restart.dump). You can then add it using self._add_file_to_current_property_instance().
self._add_property_instance_and_common_crystal_genome_keys("crystal-structure-npt",write_temp=True,write_stress=True)
self._add_file_to_current_property_instance("restart-file","restart.dump")
Note
Your Test Driver should take a fixed seed with a default value as an argument. The user can randomize this input if needed. If you run multiple simulations in your Test Driver, please derive additional seeds deterministically from the input seed.
Todo
Update kim_tools.symmetry_util.core.change_of_basis_atoms() to incorporate the functionality of kim_tools.symmetry_util.core.reduce_and_avg(), without the restriction on how the supercell is built.
Listed below are some additional considerations for writing LAMMPS MD test drivers. See https://github.com/openkim-hackathons/NPTCrystalStructure for an example Test Driver that uses LAMMPS to compute an equilibrium crystal structure under NPT conditions. We thank Philipp Hoellmer, Guanming Zhang, Tom Egg (NYU), Jason Ogbebor, and Killian Sheriff (MIT) for developing this example.
The example shows how to use an MSD check to quickly fail if a phase transition occurs, and uses the
kim-convergencepackage to ensure the simulation is equilibrated.Not all models support all LAMMPS systems of units. To work around this, use the
unit_conversion_modeoption for thekim initcommand and use the LAMMPS variables created by this command to convert to a single set of units. See https://docs.lammps.org/kim_commands.html#openkim-im-initialization-kim-init for more info.The format of a LAMMPS data file must be commensurate to the LAMMPS
atom_styledeclared. Most KIM models useatom_style atomic, however some useatom_style charge. To determine theatom_style, useself._get_supported_lammps_atom_style(). You can then pass the value returned by this function toase.io.write()to generate a LAMMPS data file with the appropriateatom_style.By default, LAMMPS will un-skew the box if it gets too tilted. This can cause
_update_nominal_parameter_values()to fail. To suppress this LAMMPS behavior, use theflip nooption withfix nptorfix deform. See https://docs.lammps.org/Howto_triclinic.html#periodicity-and-tilt-factors-for-triclinic-simulation-boxes for more info.If you are running the LAMMPS executable (e.g. using
subprocess), your Test Driver should take the name of the LAMMPS executable as an argument namedlammps_commandwith the default value"lmp". In this case, there should not be internal logic in the Test Driver for running LAMMPS with MPI -- your Driver should support being passed e.g."mpirun -np 2 lmp". This way, the user (or OpenKIM pipeline) can explicitly handle running in parallel on their machine.
Todo
LAMMPS capabilities should be better integrated into kim-tools and better examples should be provided.