Ctrl+K
Ctrl+K

Minimization in LAMMPS

Minimization in LAMMPS#

The following example will show how one performs a LAMMPS minimization calculation through AiiDA’s Python API. For the purpose of this demonstration, the optimal structure for bcc iron will is computed.

Tip

The code shown in the snippets below can be downloaded as a script, The script can be made executable and then run to execute the example calculation.

First import the required classes and functions:

from aiida.engine import run
from aiida.orm import Dict, StructureData, load_code
from aiida_lammps.data.potential import LammpsPotentialData

Then, load the code that was setup in AiiDA for lmp and get an instance of the process builder:

# Load the code configured for ``lmp``. Make sure to replace
# this string with the label used in the code setup.
code = load_code('lammps@localhost')
builder = code.get_builder()

The process builder can be used to assign and automatically validate the inputs that will be used for the calculation. One can start by defining and assigning the structure to the builder:

from ase.build import bulk
structure = StructureData(ase=bulk('Fe', 'bcc', 2.87, cubic=True))
builder.structure = structure

The crystal structure is generated by making use of the bulk method of the ASE library. The ase structure is then passed to the StructureData, which generates a node that is stored in AiiDA’s provenance graph, and then passed to the builder.

Note

The structure is constructed in such as way as to have a cubic cell with orthogonal axis as LAMMPS prefers this kind of setup.

The next step is to define the interatomic potential that will be used for this system. The interatomic potentials can be found in many repositories such as the ones from NIST and OpenKIM. In the following a potential from the OpenKIM repository will be used.

import requests
import io

# Download the potential from the repository and store it as a BytesIO object
_stream = io.BytesIO(requests.get('https://openkim.org/files/MO_546673549085_000/Fe_2.eam.fs').text.encode('ascii'))

# Set the metadata for the potential
potential_parameters = {
    'species': ['Fe'],
    'atom_style': 'atomic',
    'pair_style': 'eam/fs',
    'units': 'metal',
    'extra_tags': {
        'title': 'EAM potential (LAMMPS cubic hermite tabulation) for Fe developed by Mendelev et al. (2003) v000',
        'content_origin': 'NIST IPRP: https: // www.ctcms.nist.gov/potentials/Fe.html',
        'developer': ['Ronald E. Miller'],
        'publication_year': 2018,
    }
}

# Store the potential in an AiiDA node
potential = LammpsPotentialData.get_or_create(source=_stream,**potential_parameters)

builder.potential = potential

The potential is downloaded from the OpenKIM repository making use of the requests library and then transformed into a bytes stream via the BytesIO class from the core io library. After that one needs to define the metadata for the potential, this is necessary to be able to make sure that the potential is used only on appropriate systems and that it can be easily tracked.

Note

All the parameters in the extra_tags are not necessary to define a potential node, but they will improve the querying and tracking of the potential. For a more detailed explanation see the topic section LammpsPotentialData.

With the potential data and the metadata dictionary one can then generate a LammpsPotentialData node which can be assigned to the builder.

Then one needs to define the parameters which control how the input file for the LAMMPS calculation is generated. For a structural minimization the minimal set of parameters is the following:


# Parameters to control the input file generation
parameters = Dict({
    "control": {
        "units": "metal",
        "timestep": 1e-5
    },
    "compute":{
        "pe/atom": [{"type": [{"keyword": " ", "value": " "}], "group": "all"}],
        "ke/atom": [{"type": [{"keyword": " ", "value": " "}], "group": "all"}],
        "stress/atom": [{"type": ["NULL"], "group": "all"}],
        "pressure": [{"type": ["thermo_temp"], "group": "all"}],
    },

    "structure":{"atom_style": "atomic"},
    "thermo":{
        "printing_rate": 100,
        "thermo_printing": {
            "step": True,
            "pe": True,
            "ke": True,
            "press": True,
            "pxx": True,
            "pyy": True,
            "pzz": True,
        }
    },
    "minimize":{
        "style": "cg",
        "energy_tolerance": 1e-4,
        "force_tolerance": 1e-4,
        "max_iterations": 1000,
        "max_evaluations": 1000,
    },
})
builder.parameters = parameters

The parameters have several sections which control different behavior of the calculation:

  • control section handles global parameters for the simulation, such as the units and time step.

  • compute section specifies which parameters will be calculated and printed to file during the LAMMPS simulation (see compute command).

  • structure: controls aspects related to the structure handling in LAMMPS.

  • thermo: controls which global thermodynamic information will be calculated and written to file (see thermo command).

  • minimize: controls the structural minimization algorithm used (see minimize command).

Note

For a detailed explanation of the possible options in the parameters see the Parameters topic section.

Lastly one needs to define the computational resources needed to perform the calculation

# Run the calculation on 1 CPU and kill it if it runs longer than 1800 seconds.
# Set ``withmpi`` to ``False`` if ``pw.x`` was compiled without MPI support.
builder.metadata.options = {
    'resources': {
        'num_machines': 1,
    },
    'max_wallclock_seconds': 1800,
    'withmpi': False,
}

Now as all the needed parameters have been defined the calculation can bse launched using the process builder:

outputs, node = run.get_node(builder)

Once the calculation is finished run.get_node will return the outputs produced and the calculation node, outputs and node respectively.

The node is the entry that contains the information pertaining the calculation. It is possible to check if the calculation finished successfully (processes that return 0 are considered to be successful) by looking at its exit status:

node.exit_status

If the result is different from zero it means that a problem was encountered in the calculation. This might indicate that some output is not present, that the calculation failed due to a transitory issue, an input problem, etc.

The outputs is a dictionary containing the output nodes produced by the calculation:

print(outputs)
{
    'results': <Dict: uuid: 1dfa9349-7397-48ab-8a02-df267e3504bb (pk: 77503)>,
    'time_dependent_computes': <ArrayData: uuid: 6ec8947f-2a82-4c7a-898a-dd058d8e914e (pk: 77504)>,
    'trajectories': <LammpsTrajectory: uuid: 8755ba9d-6145-43ea-88fe-3a53a753e5eb (pk: 77505)>,
    'structure': <StructureData: uuid: 11cb6a62-01b6-464f-a737-9774f9baa9b7 (pk: 77506)>,
    'remote_folder': <RemoteData: uuid: 00f6ed4d-c026-4f1e-8d82-032b6ec4603c (pk: 77501)>,
    'retrieved': <FolderData: uuid: 7a5fea76-f007-4109-971d-082c8db38993 (pk: 77502)>
}

The results node is a dictionary that will contain the final result of the calculated thermodynamic variables as well as general information about the calculation status

print(outputs['results'].get_dict())
{
    'final_ke': 0,
    'final_pe': -8.2418066986197,
    'final_pxx': -27037.610112703,
    'final_pyy': -27037.610112703,
    'final_pzz': -27037.610112703,
    ...
}

The time_dependent_computes contains a series of numpy arrays each one representing the values the global variables as a function of time. One can get the names of the arrays by making use of the .get_arraynames() method and an individual array can be accesses by making use of .get_array('array_name') where 'array_name' is one of the strings found with the previous command.

The complete output that was written by LAMMPS to stdout, can be retrieved as follows:

results['retrieved'].base.repository.get_object_content('aiida_lammps.out')