Running MD in LAMMPS

This is a reference for how the atomate2 flows for running MD with LAMMPS can be initialized. These flows were written with solids in mind (i.e., primarily pair_style interactions with no real bond topologies), and as such are based on using the Pymatgen Structure objects as inputs.

Before running these workflows, ensure that your atomate2.yaml file in your config directory has keys LAMMPS_CMD and/or LAMMPS_MPICMD specified, and these point to a pre-compiled lammps executable.

import gzip
from importlib.resources import files as import_resource_file
from pathlib import Path

from jobflow import run_locally
from pymatgen.core import Structure
from pymatgen.io.lammps.generators import LammpsInputFile

from atomate2.lammps.flows.core import MeltQuenchThermalizeMaker
from atomate2.lammps.jobs.core import CustomLammpsMaker, LammpsNPTMaker, LammpsNVTMaker

original_force_field_file = (
    Path(import_resource_file("atomate2")).resolve()
    / ".."
    / ".."
    / "tests/test_data/lammps/Si.tersoff.gz"
)
force_field_file = (
    original_force_field_file.parent / original_force_field_file.name.split(".gz")[0]
)
with (
    gzip.open(original_force_field_file, "rb") as f_in,
    open(force_field_file, "wb") as f_out,
):
    f_out.write(f_in.read())

Running a LAMMPS simulation requires 3 files:

  1. in.lammps : input file with all the necessary parameters for the simulation

  2. forcefield.lammps : Contains all the parameters needed to construct the force field

  3. system.data : Contains data about atoms/bonding and the simulation box

In these flows, #1 can be specified by the user or taken from templates in pymatgen.io.lammps.templates. #2 is specified either as a string or a dict with the keys usually associated with forcefields, #3 is provided either as a Pymatgen Structure or a LammpsData object.

Let’s start with the forcefield file. Both representations below are equivalent:

force_field = {
    "pair_style": "tersoff",
    "pair_coeff": f"* * \
        {force_field_file} Si",
    "species": ["Si"],
}
# Can also be specified as: force_field = LammpsForceField.from_dict(force_field)
# this representation performs a basic validation of the forcefield parameters

Due to the wide vareity of forcefields and forcefield formats out there, these inputs are directly written to a forcefield.lammps file which the input script MUST include.

With the forcefield defined, we have 2 options for how to define the simulation:

  1. Use the implemented makers (NVT/NPT/NVE/Minimization/MeltQuench)

  2. Use a custom input file

Approach #1 is of use if you want to run a standard MD simulation with predefined inputs that match MD settings from other atomate2 workflows (such as ASEMD or VASPMD for example). Meanwhile, #2 is the way to go if you have a more complex MD simulation, for which you already have a pre-written input file.

Let’s first see how to do #1:

Using predefined Makers

maker = LammpsNPTMaker(force_field=force_field)

To get more control over the parameters, you must define the InputSet as follows:

settings = {
    "barostat": "berendsen",
    "start_temp": 1000,
    "end_temp": 300,
    "start_pressure": 1.0,
    "end_pressure": 1.0,
    "timestep": 0.001,
    "nsteps": 1000,
    "friction": 0.1,
}

maker.input_set_generator.update_settings(settings)

This way, all settings that have to be updated are passed in through the settings attribute. Based on the Maker/InputSetGenerator, sensible defaults are provided to avoid having to specify everything. All the default settings are present in the pymatgen.io.lammps.generators._BASE_LAMMPS_SETTINGS object.

Another thing to note: the force field can be specified as input either to the maker or the input set generator to allow for flexibility. Also: All units by default are in “metal” to better match the other solid-state MD sets.

# example structure
si_structure = Structure(
    lattice=[[0, 0, 2.73], [2.73, 0, 0], [0, 2.73, 0]],
    species=["Si", "Si"],
    coords=[[0, 0, 0], [0.5, 0.5, 0.5]],
)

Once everything is defined, make the job. If your simulation relies on additional data (such as extra force field files or bond topologies), provide that as an arguement to the .make() method of the maker. This file will be written as “extra.data” in the run directory and an additional line “include extra.data” will be written in the input file before any fixes are applied.

job = maker.make(si_structure)

Then, run the job using your prefered workflow manager. For completeness, here’s what running it locally gives:

job_output = run_locally(job, create_folders=True)
output = job_output[job.uuid][1].output

The TaskDoc by default saves a copy of the inputs and the log files generated in the simulation. It also stores the dumpfiles in the store configured with jobflow. These dumpfiles can be parsed when constructing the taskdoc and accessed via output.trajectories; but this is not done by default since lammps dump files can be excessively large. You can turn this on when defining the maker, but be warned about the runtime!

The parsed log file can be accessed as:

output.thermo_log

The full raw log file (as a string) is also stored in the TaskDoc (in the JobStore) for things that aren’t captured by the parser.

output.raw_log_file

All .dump files are parsed and stored in the JobStore by default, and can be accessed as:

output.dump_files

Custom jobs

If you instead have a custom input file, you can use the CustomLammpsMaker to run your simulation. If you want to have custom arguements in there, define them with $variables and specify the variables as dictionary.

input_file = "/path/to/input_file"
# or
input_file = "string representation of input file"
# or
input_file = LammpsInputFile.from_file(input_file)


settings = {"variable": "value"}
maker = CustomLammpsMaker(
    inputfile=input_file, force_field=force_field, settings=settings
)

A more concrete example is:

input_file = (
    "units $units \n atom_style $atom_style \n dimension 3 \n"
    "boundary p p p \n read_data input.data \n pair_style $pair_style \n"
    "include forcefield.lammps \n min_style cg \n"
    "minimize 0.0001 0.0001 1000 10000000 \n "
    "write_data run.data"
)

This inputfile does a very simple geometry relaxation. Inputs such as “units” (or anything else with a “$”) can be specified as:

settings = {"units": "metal", "pair_style": "lj/cut 2.5", "atom_style": "full"}

maker = CustomLammpsMaker(
    inputfile=input_file,
    force_field=force_field,
    settings=settings,
    validate_params=False,
    include_defaults=False,
)

job = maker.make(si_structure)
job_output = run_locally(job, create_folders=True)

This way, even custom lammps jobs can be run in high-throughput on an HPC, with access to all the benefits of atomate2 and jobflow.

If you face issues of the inputs you specify in the settings dict being incorrectly validated, manually provide validate_params=False and use_defaults=False when initializing the CustomLammpsMaker.

Melt-Quench Flow

A common usecase in MD is the generation of input geometries for a longer MD production simulation. For this task, use the MeltQuenchThermalizeMaker to automatically generate a workflow for this task.

npt_maker = LammpsNPTMaker(force_field=force_field)
nvt_maker = LammpsNVTMaker(force_field=force_field)
maker = MeltQuenchThermalizeMaker.from_temperature_steps(
    npt_maker,
    nvt_maker,
    melt_temperature=5000,
    n_steps_melt=5000,
    quench_temperature=500,
    n_steps_quench=5000,
    n_steps_thermalize=12500,
)

job = maker.make(si_structure)
job_output = run_locally(job, create_folders=True)