Showing 35 of 118 files from the diff.
Other files ignored by Codecov
.gitignore has changed.
setup.py has changed.
pytest.ini was deleted.
.travis.yml was deleted.
MANIFEST.in has changed.
.flake8 was deleted.
setup.json has changed.
pyproject.toml has changed.
README.md has changed.
tox.ini is new.

@@ -1,15 +1,12 @@
Loading
1 1
import pytest
2 -
import six
2 +
3 3
from aiida_lammps.common.generate_structure import generate_lammps_structure
4 4
5 5
6 -
@pytest.mark.parametrize('structure', [
7 -
    "Fe",
8 -
    "pyrite",
9 -
    "fes_cubic-zincblende",
10 -
    "greigite"
11 -
])
6 +
@pytest.mark.parametrize(
7 +
    "structure", ["Fe", "pyrite", "fes_cubic-zincblende", "greigite"]
8 +
)
12 9
def test_generate(db_test_app, get_structure_data, structure, file_regression):
13 10
    structure = get_structure_data(structure)
14 11
    text, transform = generate_lammps_structure(structure, round_dp=8)
15 -
    file_regression.check(six.ensure_text(text))
12 +
    file_regression.check(text)

@@ -1,68 +1,75 @@
Loading
1 -
from aiida.orm import ArrayData, Dict, StructureData
1 +
import traceback
2 2
3 +
from aiida.orm import Dict
4 +
5 +
from aiida_lammps.common.raw_parsers import get_units_dict
6 +
from aiida_lammps.data.trajectory import LammpsTrajectory
3 7
from aiida_lammps.parsers.lammps.base import LAMMPSBaseParser
4 -
from aiida_lammps.common.raw_parsers import read_lammps_positions_and_forces_txt, get_units_dict
5 8
6 9
7 10
class OptimizeParser(LAMMPSBaseParser):
8 -
    """
9 -
    Simple Optimize Parser for LAMMPS.
10 -
    """
11 +
    """Parser for LAMMPS optimization calculation."""
11 12
12 13
    def __init__(self, node):
13 -
        """
14 -
        Initialize the instance of Optimize LammpsParser
15 -
        """
14 +
        """Initialize the instance of Optimize Lammps Parser."""
16 15
        super(OptimizeParser, self).__init__(node)
17 16
18 17
    def parse(self, **kwargs):
19 -
        """
20 -
        Parses the datafolder, stores results.
21 -
        """
22 -
        resources, exit_code = self.get_parsing_resources(kwargs)
23 -
        if exit_code is not None:
24 -
            return exit_code
25 -
        trajectory_filename, trajectory_filepath, info_filepath = resources
18 +
        """Parses the datafolder, stores results."""
19 +
        resources = self.get_parsing_resources(kwargs, traj_in_temp=True)
20 +
        if resources.exit_code is not None:
21 +
            return resources.exit_code
26 22
27 -
        log_data, exit_code = self.parse_log_file(compute_stress=True)
23 +
        log_data, exit_code = self.parse_log_file()
28 24
        if exit_code is not None:
29 25
            return exit_code
30 26
31 -
        trajectory_txt = self.retrieved.get_object_content(trajectory_filename)
32 -
        if not trajectory_txt:
33 -
            self.logger.error("trajectory file empty")
34 -
            return self.exit_codes.ERROR_TRAJ_PARSING
35 -
        positions, forces, charges, symbols, cell2 = read_lammps_positions_and_forces_txt(
36 -
            trajectory_txt)
37 -
38 -
        # save optimized structure into node
39 -
        structure = StructureData(cell=log_data["cell"])
40 -
        for i, position in enumerate(positions[-1]):
41 -
            structure.append_atom(position=position.tolist(),
42 -
                                  symbols=symbols[i])
43 -
        self.out('structure', structure)
44 -
45 -
        # save forces and stresses into node
46 -
        array_data = ArrayData()
47 -
        array_data.set_array('forces', forces)
48 -
        array_data.set_array('stress', log_data["stress"])
49 -
        array_data.set_array('positions', positions)
50 -
        if charges is not None:
51 -
            array_data.set_array('charges', charges)
52 -
        self.out('arrays', array_data)
27 +
        traj_error = None
28 +
        if not resources.traj_paths:
29 +
            traj_error = self.exit_codes.ERROR_TRAJ_FILE_MISSING
30 +
        else:
31 +
            try:
32 +
                trajectory_data = LammpsTrajectory(
33 +
                    resources.traj_paths[0],
34 +
                    aliases={
35 +
                        "stresses": ["c_stpa[{}]".format(i + 1) for i in range(6)],
36 +
                        "forces": ["fx", "fy", "fz"],
37 +
                    },
38 +
                )
39 +
                self.out("trajectory_data", trajectory_data)
40 +
                self.out(
41 +
                    "structure",
42 +
                    trajectory_data.get_step_structure(
43 +
                        -1, original_structure=self.node.inputs.structure
44 +
                    ),
45 +
                )
46 +
            except Exception as err:
47 +
                traceback.print_exc()
48 +
                self.logger.error(str(err))
49 +
                traj_error = self.exit_codes.ERROR_TRAJ_PARSING
53 50
54 51
        # save results into node
55 52
        output_data = log_data["data"]
56 -
        if 'units_style' in output_data:
57 -
            output_data.update(get_units_dict(output_data['units_style'],
58 -
                                              ["energy", "force", "distance", "pressure"]))
53 +
        if "units_style" in output_data:
54 +
            output_data.update(
55 +
                get_units_dict(
56 +
                    output_data["units_style"],
57 +
                    ["energy", "force", "distance", "pressure"],
58 +
                )
59 +
            )
59 60
            output_data["stress_units"] = output_data.pop("pressure_units")
60 61
        else:
61 62
            self.logger.warning("units missing in log")
62 63
        self.add_warnings_and_errors(output_data)
63 64
        self.add_standard_info(output_data)
64 65
        parameters_data = Dict(dict=output_data)
65 -
        self.out('results', parameters_data)
66 +
        self.out("results", parameters_data)
66 67
67 68
        if output_data["errors"]:
68 69
            return self.exit_codes.ERROR_LAMMPS_RUN
70 +
71 +
        if traj_error:
72 +
            return traj_error
73 +
74 +
        if not log_data.get("found_end", False):
75 +
            return self.exit_codes.ERROR_RUN_INCOMPLETE

@@ -0,0 +1,117 @@
Loading
1 +
from collections import namedtuple
2 +
3 +
from aiida.orm import StructureData
4 +
import numpy as np
5 +
6 +
TRAJ_BLOCK = namedtuple(
7 +
    "TRAJ_BLOCK", ["lines", "timestep", "natoms", "cell", "pbc", "atom_fields"]
8 +
)
9 +
10 +
11 +
def iter_step_lines(file_obj):
12 +
    step_content = None
13 +
    init_line = 0
14 +
    for i, line in enumerate(file_obj):
15 +
        if "ITEM: TIMESTEP" in line:
16 +
            if step_content is not None:
17 +
                yield init_line, step_content
18 +
            init_line = i + 1
19 +
            step_content = []
20 +
        if step_content is not None:
21 +
            step_content.append(line.strip())
22 +
    if step_content is not None:
23 +
        yield init_line, step_content
24 +
25 +
26 +
def parse_step(lines, intial_line=0):
27 +
    if "ITEM: TIMESTEP" not in lines[0]:
28 +
        raise IOError("expected line {} to be TIMESTEP".format(intial_line))
29 +
    if "ITEM: NUMBER OF ATOMS" not in lines[2]:
30 +
        raise IOError("expected line {} to be NUMBER OF ATOMS".format(intial_line + 2))
31 +
    if "ITEM: BOX BOUNDS xy xz yz" not in lines[4]:
32 +
        raise IOError(
33 +
            "expected line {} to be BOX BOUNDS xy xz yz".format(intial_line + 4)
34 +
        )
35 +
        # TODO handle case when xy xz yz not present -> orthogonal box
36 +
    if "ITEM: ATOMS" not in lines[8]:
37 +
        raise IOError("expected line {} to be ATOMS".format(intial_line + 8))
38 +
    timestep = int(lines[1])
39 +
    number_of_atoms = int(lines[3])
40 +
41 +
    # each pbc contains two letters <lo><hi> such that:
42 +
    # p = periodic, f = fixed, s = shrink wrap, m = shrink wrapped with a minimum value
43 +
    pbc = lines[4].split()[6:]
44 +
45 +
    bounds = [line.split() for line in lines[5:8]]
46 +
    bounds = np.array(bounds, dtype=float)
47 +
    if bounds.shape[1] == 2:
48 +
        bounds = np.append(bounds, np.array([0, 0, 0])[None].T, axis=1)
49 +
50 +
    xy = bounds[0, 2]
51 +
    xz = bounds[1, 2]
52 +
    yz = bounds[2, 2]
53 +
54 +
    xlo = bounds[0, 0] - np.min([0.0, xy, xz, xy + xz])
55 +
    xhi = bounds[0, 1] - np.max([0.0, xy, xz, xy + xz])
56 +
    ylo = bounds[1, 0] - np.min([0.0, yz])
57 +
    yhi = bounds[1, 1] - np.max([0.0, yz])
58 +
    zlo = bounds[2, 0]
59 +
    zhi = bounds[2, 1]
60 +
61 +
    super_cell = np.array([[xhi - xlo, xy, xz], [0, yhi - ylo, yz], [0, 0, zhi - zlo]])
62 +
    cell = super_cell.T
63 +
    field_names = lines[8].split()[2:]
64 +
    fields = []
65 +
    for i in range(number_of_atoms):
66 +
        fields.append(lines[9 + i].split())
67 +
    atom_fields = {n: v.tolist() for n, v in zip(field_names, np.array(fields).T)}
68 +
69 +
    return TRAJ_BLOCK(lines, timestep, number_of_atoms, cell, pbc, atom_fields)
70 +
71 +
72 +
def iter_trajectories(file_obj):
73 +
    """Parse a LAMMPS Trajectory file, yielding data for each time step."""
74 +
    for line_num, lines in iter_step_lines(file_obj):
75 +
        yield parse_step(lines, line_num)
76 +
77 +
78 +
def create_structure(
79 +
    traj_block,
80 +
    symbol_field="element",
81 +
    position_fields=("x", "y", "z"),
82 +
    original_structure=None,
83 +
):
84 +
    symbols = traj_block.atom_fields[symbol_field]
85 +
    positions = np.array(
86 +
        [traj_block.atom_fields[f] for f in position_fields], dtype=float
87 +
    ).T
88 +
89 +
    if original_structure is not None:
90 +
        kind_names = original_structure.get_site_kindnames()
91 +
        kind_symbols = [original_structure.get_kind(n).symbol for n in kind_names]
92 +
        if symbols != kind_symbols:
93 +
            raise ValueError(
94 +
                "original_structure has different symbols:: {} != {}".format(
95 +
                    kind_symbols, symbols
96 +
                )
97 +
            )
98 +
        structure = original_structure.clone()
99 +
        structure.reset_cell(traj_block.cell)
100 +
        structure.reset_sites_positions(positions)
101 +
102 +
        return structure
103 +
104 +
    pbcs = []
105 +
    for pbc in traj_block.pbc:
106 +
        if pbc == "pp":
107 +
            pbcs.append(True)
108 +
        elif pbc == "ff":
109 +
            pbcs.append(False)
110 +
        else:
111 +
            raise NotImplementedError("pbc = {}".format(traj_block.pbc))
112 +
113 +
    structure = StructureData(cell=traj_block.cell, pbc=pbcs)
114 +
    for symbol, position in zip(symbols, positions):
115 +
        structure.append_atom(position=position, symbols=symbol)
116 +
117 +
    return structure

@@ -1,43 +1,6 @@
Loading
1 -
import json
2 -
import os
3 -
4 -
import jsonschema
5 -
6 -
7 -
def read_schema(name):
8 -
    """read and return an json schema
9 -
10 -
    :param name: <name>.schema.json
11 -
12 -
    :return: dictionary
13 -
    """
14 -
    dirpath = os.path.dirname(os.path.realpath(__file__))
15 -
    jpath = os.path.join(dirpath, "{}.schema.json".format(name))
16 -
    with open(jpath) as jfile:
17 -
        schema = json.load(jfile)
18 -
    return schema
19 -
20 -
21 -
def validate_with_json(data, name):
22 -
    """ validate json-type data against a schema
23 -
24 -
    :param data: dictionary
25 -
    :param name: <name>.schema.json
26 -
    """
27 -
    schema = read_schema(name)
28 -
    validator = jsonschema.Draft4Validator
29 -
30 -
    # by default, only validates lists
31 -
    validator(schema, types={"array": (list, tuple)}).validate(data)
32 -
33 -
34 -
def validate_with_dict(data, schema):
35 -
    """ validate json-type data against a schema
36 -
37 -
    :param data: dictionary
38 -
    :param schema: dictionary
39 -
    """
40 -
    validator = jsonschema.Draft4Validator
41 -
42 -
    # by default, only validates lists
43 -
    validator(schema, types={"array": (list, tuple)}).validate(data)
1 +
"""Package for validating JSON objects against schemas."""
2 +
from aiida_lammps.validation.utils import (  # noqa: F401
3 +
    load_schema,
4 +
    load_validator,
5 +
    validate_against_schema,
6 +
)
44 7
imilarity index 99%
45 8
ename from aiida_lammps/validation/force.schema.json
46 9
ename to aiida_lammps/validation/schemas/force.schema.json

@@ -1,14 +1,14 @@
Loading
1 -
from aiida.engine.calculation.job import CalcJob
2 -
from aiida.orm.nodes.parameter import Dict
3 -
from aiida.orm.nodes.structure import StructureData
4 -
from aiida.orm.nodes.array.trajectory import TrajectoryData
5 -
from aiida.orm.nodes.array import ArrayData
6 -
from aiida.common.exceptions import InputValidationError
7 1
from aiida.common.datastructures import CalcInfo, CodeInfo
2 +
from aiida.common.exceptions import InputValidationError
8 3
from aiida.common.utils import classproperty
4 +
from aiida.engine import CalcJob
5 +
from aiida.orm import ArrayData, StructureData, TrajectoryData
6 +
from aiida_phonopy.common.raw_parsers import get_force_constants, get_poscar_txt
9 7
10 -
from aiida_phonopy.common.raw_parsers import get_FORCE_CONSTANTS_txt, get_poscar_txt
11 -
from aiida_lammps.common.generate_input_files import get_trajectory_txt, parameters_to_input_file
8 +
from aiida_lammps.common.generate_input_files import (
9 +
    get_trajectory_txt,
10 +
    parameters_to_input_file,
11 +
)
12 12
13 13
14 14
class DynaphopyCalculation(CalcJob):
@@ -21,17 +21,16 @@
Loading
21 21
    def _init_internal_params(self):
22 22
        super(DynaphopyCalculation, self)._init_internal_params()
23 23
24 -
        self._INPUT_FILE_NAME = 'input_dynaphopy'
25 -
        self._INPUT_TRAJECTORY = 'trajectory'
26 -
        self._INPUT_CELL = 'POSCAR'
27 -
        self._INPUT_FORCE_CONSTANTS = 'FORCE_CONSTANTS'
28 -
29 -
        self._OUTPUT_FORCE_CONSTANTS = 'FORCE_CONSTANTS_OUT'
30 -
        self._OUTPUT_FILE_NAME = 'OUTPUT'
31 -
        self._OUTPUT_QUASIPARTICLES = 'quasiparticles_data.yaml'
24 +
        self._INPUT_FILE_NAME = "input_dynaphopy"
25 +
        self._INPUT_TRAJECTORY = "trajectory"
26 +
        self._INPUT_CELL = "POSCAR"
27 +
        self._INPUT_FORCE_CONSTANTS = "FORCE_CONSTANTS"
32 28
33 -
        self._default_parser = 'dynaphopy'
29 +
        self._OUTPUT_FORCE_CONSTANTS = "FORCE_CONSTANTS_OUT"
30 +
        self._OUTPUT_FILE_NAME = "OUTPUT"
31 +
        self._OUTPUT_QUASIPARTICLES = "quasiparticles_data.yaml"
34 32
33 +
        self._default_parser = "dynaphopy"
35 34
36 35
    @classproperty
37 36
    def _use_methods(cls):
@@ -39,38 +38,46 @@
Loading
39 38
        Additional use_* methods for the namelists class.
40 39
        """
41 40
        retdict = JobCalculation._use_methods
42 -
        retdict.update({
43 -
            "parameters": {
44 -
               'valid_types': ParameterData,
45 -
               'additional_parameter': None,
46 -
               'linkname': 'parameters',
47 -
               'docstring': ("Use a node that specifies the dynaphopy input "
48 -
                             "for the namelists"),
49 -
               },
50 -
            "trajectory": {
51 -
               'valid_types': TrajectoryData,
52 -
               'additional_parameter': None,
53 -
               'linkname': 'trajectory',
54 -
               'docstring': ("Use a node that specifies the trajectory data "
55 -
                             "for the namelists"),
56 -
               },
57 -
            "force_constants": {
58 -
               'valid_types': ArrayData,
59 -
               'additional_parameter': None,
60 -
               'linkname': 'force_constants',
61 -
               'docstring': ("Use a node that specifies the force_constants "
62 -
                             "for the namelists"),
63 -
               },
64 -
            "structure": {
65 -
               'valid_types': StructureData,
66 -
               'additional_parameter': None,
67 -
               'linkname': 'structure',
68 -
               'docstring': "Use a node for the structure",
69 -
               },
70 -
         })
41 +
        retdict.update(
42 +
            {
43 +
                "parameters": {
44 +
                    "valid_types": ParameterData,
45 +
                    "additional_parameter": None,
46 +
                    "linkname": "parameters",
47 +
                    "docstring": (
48 +
                        "Use a node that specifies the dynaphopy input "
49 +
                        "for the namelists"
50 +
                    ),
51 +
                },
52 +
                "trajectory": {
53 +
                    "valid_types": TrajectoryData,
54 +
                    "additional_parameter": None,
55 +
                    "linkname": "trajectory",
56 +
                    "docstring": (
57 +
                        "Use a node that specifies the trajectory data "
58 +
                        "for the namelists"
59 +
                    ),
60 +
                },
61 +
                "force_constants": {
62 +
                    "valid_types": ArrayData,
63 +
                    "additional_parameter": None,
64 +
                    "linkname": "force_constants",
65 +
                    "docstring": (
66 +
                        "Use a node that specifies the force_constants "
67 +
                        "for the namelists"
68 +
                    ),
69 +
                },
70 +
                "structure": {
71 +
                    "valid_types": StructureData,
72 +
                    "additional_parameter": None,
73 +
                    "linkname": "structure",
74 +
                    "docstring": "Use a node for the structure",
75 +
                },
76 +
            }
77 +
        )
71 78
        return retdict
72 79
73 -
    def _prepare_for_submission(self,tempfolder, inputdict):
80 +
    def _prepare_for_submission(self, tempfolder, inputdict):
74 81
        """
75 82
        This is the routine to be called when you want to create
76 83
        the input files and related stuff with a plugin.
@@ -82,37 +89,37 @@
Loading
82 89
        """
83 90
84 91
        try:
85 -
            parameters_data = inputdict.pop(self.get_linkname('parameters'))
92 +
            parameters_data = inputdict.pop(self.get_linkname("parameters"))
86 93
        except KeyError:
87 94
            pass
88 -
            #raise InputValidationError("No parameters specified for this "
95 +
            # raise InputValidationError("No parameters specified for this "
89 96
            #                           "calculation")
90 97
        if not isinstance(parameters_data, ParameterData):
91 -
            raise InputValidationError("parameters is not of type "
92 -
                                       "ParameterData")
98 +
            raise InputValidationError("parameters is not of type " "ParameterData")
93 99
94 100
        try:
95 -
            structure = inputdict.pop(self.get_linkname('structure'))
101 +
            structure = inputdict.pop(self.get_linkname("structure"))
96 102
        except KeyError:
97 103
            raise InputValidationError("no structure is specified for this calculation")
98 104
99 105
        try:
100 -
            trajectory = inputdict.pop(self.get_linkname('trajectory'))
106 +
            trajectory = inputdict.pop(self.get_linkname("trajectory"))
101 107
        except KeyError:
102 108
            raise InputValidationError("trajectory is specified for this calculation")
103 109
104 110
        try:
105 -
            force_constants = inputdict.pop(self.get_linkname('force_constants'))
111 +
            force_constants = inputdict.pop(self.get_linkname("force_constants"))
106 112
        except KeyError:
107 -
            raise InputValidationError("no force_constants is specified for this calculation")
113 +
            raise InputValidationError(
114 +
                "no force_constants is specified for this calculation"
115 +
            )
108 116
109 117
        try:
110 -
            code = inputdict.pop(self.get_linkname('code'))
118 +
            code = inputdict.pop(self.get_linkname("code"))
111 119
        except KeyError:
112 120
            raise InputValidationError("no code is specified for this calculation")
113 121
114 -
115 -
        time_step = trajectory.get_times()[1]-trajectory.get_times()[0]
122 +
        time_step = trajectory.get_times()[1] - trajectory.get_times()[0]
116 123
117 124
        ##############################
118 125
        # END OF INITIAL INPUT CHECK #
@@ -122,32 +129,32 @@
Loading
122 129
123 130
        cell_txt = get_poscar_txt(structure)
124 131
        input_txt = parameters_to_input_file(parameters_data)
125 -
        force_constants_txt = get_FORCE_CONSTANTS_txt(force_constants)
132 +
        force_constants_txt = get_force_constants(force_constants)
126 133
        trajectory_txt = get_trajectory_txt(trajectory)
127 134
128 135
        # =========================== dump to file =============================
129 136
130 137
        input_filename = tempfolder.get_abs_path(self._INPUT_FILE_NAME)
131 -
        with open(input_filename, 'w') as infile:
138 +
        with open(input_filename, "w") as infile:
132 139
            infile.write(input_txt)
133 140
134 141
        cell_filename = tempfolder.get_abs_path(self._INPUT_CELL)
135 -
        with open(cell_filename, 'w') as infile:
142 +
        with open(cell_filename, "w") as infile:
136 143
            infile.write(cell_txt)
137 144
138 145
        force_constants_filename = tempfolder.get_abs_path(self._INPUT_FORCE_CONSTANTS)
139 -
        with open(force_constants_filename, 'w') as infile:
146 +
        with open(force_constants_filename, "w") as infile:
140 147
            infile.write(force_constants_txt)
141 148
142 149
        trajectory_filename = tempfolder.get_abs_path(self._INPUT_TRAJECTORY)
143 -
        with open(trajectory_filename, 'w') as infile:
150 +
        with open(trajectory_filename, "w") as infile:
144 151
            infile.write(trajectory_txt)
145 152
146 153
        # ============================ calcinfo ================================
147 154
148 155
        local_copy_list = []
149 156
        remote_copy_list = []
150 -
    #    additional_retrieve_list = settings_dict.pop("ADDITIONAL_RETRIEVE_LIST",[])
157 +
        #    additional_retrieve_list = settings_dict.pop("ADDITIONAL_RETRIEVE_LIST",[])
151 158
152 159
        calcinfo = CalcInfo()
153 160
@@ -157,23 +164,37 @@
Loading
157 164
        calcinfo.remote_copy_list = remote_copy_list
158 165
159 166
        # Retrieve files
160 -
        calcinfo.retrieve_list = [self._OUTPUT_FILE_NAME,
161 -
                                  self._OUTPUT_FORCE_CONSTANTS,
162 -
                                  self._OUTPUT_QUASIPARTICLES]
167 +
        calcinfo.retrieve_list = [
168 +
            self._OUTPUT_FILE_NAME,
169 +
            self._OUTPUT_FORCE_CONSTANTS,
170 +
            self._OUTPUT_QUASIPARTICLES,
171 +
        ]
163 172
164 173
        codeinfo = CodeInfo()
165 -
        codeinfo.cmdline_params = [self._INPUT_FILE_NAME, self._INPUT_TRAJECTORY,
166 -
                                   '-ts', '{}'.format(time_step), '--silent',
167 -
                                   '-sfc', self._OUTPUT_FORCE_CONSTANTS, '-thm', # '--resolution 0.01',
168 -
                                   '-psm','2', '--normalize_dos', '-sdata']
169 -
170 -
        if 'temperature' in parameters_data.get_dict():
171 -
            codeinfo.cmdline_params.append('--temperature')
172 -
            codeinfo.cmdline_params.append('{}'.format(parameters_data.dict.temperature))
173 -
174 -
        if 'md_commensurate' in parameters_data.get_dict():
174 +
        codeinfo.cmdline_params = [
175 +
            self._INPUT_FILE_NAME,
176 +
            self._INPUT_TRAJECTORY,
177 +
            "-ts",
178 +
            "{}".format(time_step),
179 +
            "--silent",
180 +
            "-sfc",
181 +
            self._OUTPUT_FORCE_CONSTANTS,
182 +
            "-thm",  # '--resolution 0.01',
183 +
            "-psm",
184 +
            "2",
185 +
            "--normalize_dos",
186 +
            "-sdata",
187 +
        ]
188 +
189 +
        if "temperature" in parameters_data.get_dict():
190 +
            codeinfo.cmdline_params.append("--temperature")
191 +
            codeinfo.cmdline_params.append(
192 +
                "{}".format(parameters_data.dict.temperature)
193 +
            )
194 +
195 +
        if "md_commensurate" in parameters_data.get_dict():
175 196
            if parameters_data.dict.md_commensurate:
176 -
                codeinfo.cmdline_params.append('--MD_commensurate')
197 +
                codeinfo.cmdline_params.append("--MD_commensurate")
177 198
178 199
        codeinfo.stdout_name = self._OUTPUT_FILE_NAME
179 200
        codeinfo.code_uuid = code.uuid

@@ -1,29 +1,175 @@
Loading
1 +
"""Module for parsing REAXFF input files.
2 +
3 +
Note: this module is copied directly from aiida-crystal17 v0.10.0b5
4 +
"""
1 5
import copy
2 6
import re
3 -
from decimal import Decimal
4 -
5 -
import numpy as np
6 -
7 -
8 -
# TODO can X be in middle of species?
9 -
10 -
11 -
def skip(f, numlines):
12 -
    for _ in range(numlines):
13 -
        f.readline()
14 -
15 -
16 -
def readnumline(f, vtype):
17 -
    line = f.readline().split()
18 -
    if line[1] != '!':
19 -
        raise Exception(
20 -
            'not on a line containing ! (as expected) while reading {0}'.
21 -
            format(vtype))
22 -
    return int(line[0])
23 7
8 +
from aiida_lammps.validation import validate_against_schema
9 +
10 +
INDEX_SEP = "-"
11 +
12 +
KEYS_GLOBAL = (
13 +
    "reaxff0_boc1",
14 +
    "reaxff0_boc2",
15 +
    "reaxff3_coa2",
16 +
    "Triple bond stabilisation 1",
17 +
    "Triple bond stabilisation 2",
18 +
    "C2-correction",
19 +
    "reaxff0_ovun6",
20 +
    "Triple bond stabilisation",
21 +
    "reaxff0_ovun7",
22 +
    "reaxff0_ovun8",
23 +
    "Triple bond stabilization energy",
24 +
    "Lower Taper-radius",
25 +
    "Upper Taper-radius",
26 +
    "reaxff2_pen2",
27 +
    "reaxff0_val7",
28 +
    "reaxff0_lp1",
29 +
    "reaxff0_val9",
30 +
    "reaxff0_val10",
31 +
    "Not used 2",
32 +
    "reaxff0_pen2",
33 +
    "reaxff0_pen3",
34 +
    "reaxff0_pen4",
35 +
    "Not used 3",
36 +
    "reaxff0_tor2",
37 +
    "reaxff0_tor3",
38 +
    "reaxff0_tor4",
39 +
    "Not used 4",
40 +
    "reaxff0_cot2",
41 +
    "reaxff0_vdw1",
42 +
    "bond order cutoff",
43 +
    "reaxff3_coa4",
44 +
    "reaxff0_ovun4",
45 +
    "reaxff0_ovun3",
46 +
    "reaxff0_val8",
47 +
    "Not used 5",
48 +
    "Not used 6",
49 +
    "Not used 7",
50 +
    "Not used 8",
51 +
    "reaxff3_coa3",
52 +
)
53 +
54 +
# TODO some variables lammps sets as global are actually species dependant in GULP, how to handle these?
55 +
56 +
KEYS_1BODY = (
57 +
    "reaxff1_radii1",
58 +
    "reaxff1_valence1",
59 +
    "mass",
60 +
    "reaxff1_morse3",
61 +
    "reaxff1_morse2",
62 +
    "reaxff_gamma",
63 +
    "reaxff1_radii2",
64 +
    "reaxff1_valence3",
65 +
    "reaxff1_morse1",
66 +
    "reaxff1_morse4",
67 +
    "reaxff1_valence4",
68 +
    "reaxff1_under",
69 +
    "dummy1",
70 +
    "reaxff_chi",
71 +
    "reaxff_mu",
72 +
    "dummy2",
73 +
    "reaxff1_radii3",
74 +
    "reaxff1_lonepair2",
75 +
    "dummy3",
76 +
    "reaxff1_over2",
77 +
    "reaxff1_over1",
78 +
    "reaxff1_over3",
79 +
    "dummy4",
80 +
    "dummy5",
81 +
    "reaxff1_over4",
82 +
    "reaxff1_angle1",
83 +
    "dummy11",
84 +
    "reaxff1_valence2",
85 +
    "reaxff1_angle2",
86 +
    "dummy6",
87 +
    "dummy7",
88 +
    "dummy8",
89 +
)
90 +
91 +
KEYS_2BODY_BONDS = (
92 +
    "reaxff2_bond1",
93 +
    "reaxff2_bond2",
94 +
    "reaxff2_bond3",
95 +
    "reaxff2_bond4",
96 +
    "reaxff2_bo5",
97 +
    "reaxff2_bo7",
98 +
    "reaxff2_bo6",
99 +
    "reaxff2_over",
100 +
    "reaxff2_bond5",
101 +
    "reaxff2_bo3",
102 +
    "reaxff2_bo4",
103 +
    "dummy1",
104 +
    "reaxff2_bo1",
105 +
    "reaxff2_bo2",
106 +
    "reaxff2_bo8",
107 +
    "reaxff2_pen1",
108 +
)
109 +
110 +
KEYS_2BODY_OFFDIAG = [
111 +
    "reaxff2_morse1",
112 +
    "reaxff2_morse3",
113 +
    "reaxff2_morse2",
114 +
    "reaxff2_morse4",
115 +
    "reaxff2_morse5",
116 +
    "reaxff2_morse6",
117 +
]
24 118
25 -
def split_numbers(string, as_decimal=False):
26 -
    """ get a list of numbers from a string (even with no spacing)
119 +
KEYS_3BODY_ANGLES = (
120 +
    "reaxff3_angle1",
121 +
    "reaxff3_angle2",
122 +
    "reaxff3_angle3",
123 +
    "reaxff3_coa1",
124 +
    "reaxff3_angle5",
125 +
    "reaxff3_penalty",
126 +
    "reaxff3_angle4",
127 +
)
128 +
129 +
KEYS_3BODY_HBOND = (
130 +
    "reaxff3_hbond1",
131 +
    "reaxff3_hbond2",
132 +
    "reaxff3_hbond3",
133 +
    "reaxff3_hbond4",
134 +
)
135 +
136 +
KEYS_4BODY_TORSION = (
137 +
    "reaxff4_torsion1",
138 +
    "reaxff4_torsion2",
139 +
    "reaxff4_torsion3",
140 +
    "reaxff4_torsion4",
141 +
    "reaxff4_torsion5",
142 +
    "dummy1",
143 +
    "dummy2",
144 +
)
145 +
146 +
DEFAULT_TOLERANCES = (
147 +
    # ReaxFF angle/torsion bond order threshold,
148 +
    # for bond orders in valence, penalty and 3-body conjugation
149 +
    # GULP default: 0.001
150 +
    ("anglemin", 0.001),
151 +
    # ReaxFF bond order double product threshold,
152 +
    # for the product of bond orders (1-2 x 2-3, where 2 = pivot)
153 +
    # Hard coded to 0.001 in original code, but this leads to discontinuities
154 +
    # GULP default: 0.000001
155 +
    ("angleprod", 0.00001),
156 +
    # ReaxFF hydrogen-bond bond order threshold
157 +
    # Hard coded to 0.01 in original code.
158 +
    # GULP default: 0.01
159 +
    ("hbondmin", 0.01),
160 +
    # ReaxFF H-bond cutoff
161 +
    # Hard coded to 7.5 Ang in original code.
162 +
    # GULP default: 7.5
163 +
    ("hbonddist", 7.5),
164 +
    # ReaxFF bond order triple product threshold,
165 +
    # for the product of bond orders (1-2 x 2-3 x 3-4)
166 +
    # GULP default: 0.000000001
167 +
    ("torsionprod", 0.00001),
168 +
)
169 +
170 +
171 +
def split_numbers(string):
172 +
    """Get a list of numbers from a string (even with no spacing).
27 173
28 174
    :type string: str
29 175
    :type as_decimal: bool
@@ -52,673 +198,463 @@
Loading
52 198
    [0.001, -2.0]
53 199
54 200
    """
55 -
    _match_number = re.compile(
56 -
        '-?\\ *[0-9]+\\.?[0-9]*(?:[Ee]\\ *[+-]?\\ *[0-9]+)?')
201 +
    _match_number = re.compile("-?\\ *[0-9]+\\.?[0-9]*(?:[Ee]\\ *[+-]?\\ *[0-9]+)?")
57 202
    string = string.replace(" .", " 0.")
58 203
    string = string.replace("-.", "-0.")
59 -
    return [
60 -
        Decimal(s) if as_decimal else float(s)
61 -
        for s in re.findall(_match_number, string)
62 -
    ]
63 -
64 -
65 -
def readlist(f):
66 -
    return split_numbers(f.readline())
67 -
68 -
69 -
def getvals(key, vals, names):
70 -
    out = []
71 -
    for name in names:
72 -
        out.append(vals[key.index(name)])
73 -
    return out
74 -
75 -
76 -
def transpose(lst):
77 -
    return list(map(list, zip(*lst)))
78 -
79 -
80 -
_paramkeys = [
81 -
    'Overcoordination 1', 'Overcoordination 2', 'Valency angle conjugation 1',
82 -
    'Triple bond stabilisation 1', 'Triple bond stabilisation 2',
83 -
    'C2-correction', 'Undercoordination 1', 'Triple bond stabilisation',
84 -
    'Undercoordination 2', 'Undercoordination 3',
85 -
    'Triple bond stabilization energy', 'Lower Taper-radius',
86 -
    'Upper Taper-radius', 'Not used 1', 'Valency undercoordination',
87 -
    'Valency angle/lone pair', 'Valency angle 1', 'Valency angle 2',
88 -
    'Not used 2', 'Double bond/angle', 'Double bond/angle: overcoord 1',
89 -
    'Double bond/angle: overcoord 2', 'Not used 3', 'Torsion/BO',
90 -
    'Torsion overcoordination 1', 'Torsion overcoordination 2', 'Not used 4',
91 -
    'Conjugation', 'vdWaals shielding', 'bond order cutoff',
92 -
    'Valency angle conjugation 2', 'Valency overcoordination 1',
93 -
    'Valency overcoordination 2', 'Valency/lone pair', 'Not used 5',
94 -
    'Not used 6', 'Not used 7', 'Not used 8', 'Valency angle conjugation 3'
95 -
]
96 -
97 -
_speckeys = [
98 -
    'idx', 'symbol', 'reaxff1_radii1', 'reaxff1_valence1', 'mass',
99 -
    'reaxff1_morse3', 'reaxff1_morse2', 'reaxff_gamma', 'reaxff1_radii2',
100 -
    'reaxff1_valence3', 'reaxff1_morse1', 'reaxff1_morse4', 'reaxff1_valence4',
101 -
    'reaxff1_under', 'dummy1', 'reaxff_chi', 'reaxff_mu', 'dummy2',
102 -
    'reaxff1_radii3', 'reaxff1_lonepair2', 'dummy3', 'reaxff1_over2',
103 -
    'reaxff1_over1', 'reaxff1_over3', 'dummy4', 'dummy5', 'reaxff1_over4',
104 -
    'reaxff1_angle1', 'dummy11', 'reaxff1_valence2', 'reaxff1_angle2',
105 -
    'dummy6', 'dummy7', 'dummy8'
106 -
]
204 +
    return [float(s) for s in re.findall(_match_number, string)]
107 205
108 -
_bondkeys = [
109 -
    'idx1', 'idx2', 'reaxff2_bond1', 'reaxff2_bond2', 'reaxff2_bond3',
110 -
    'reaxff2_bond4', 'reaxff2_bo5', 'reaxff2_bo7', 'reaxff2_bo6',
111 -
    'reaxff2_over', 'reaxff2_bond5', 'reaxff2_bo3', 'reaxff2_bo4', 'dummy1',
112 -
    'reaxff2_bo1', 'reaxff2_bo2', 'reaxff2_bo8', 'reaxff2_bo9'
113 -
]
114 -
115 -
_odkeys = [
116 -
    'idx1', 'idx2', 'reaxff2_morse1', 'reaxff2_morse3', 'reaxff2_morse2',
117 -
    'reaxff2_morse4', 'reaxff2_morse5', 'reaxff2_morse6'
118 -
]
119 206
120 -
_anglekeys = [
121 -
    'idx1', 'idx2', 'idx3', 'reaxff3_angle1', 'reaxff3_angle2',
122 -
    'reaxff3_angle3', 'reaxff3_conj', 'reaxff3_angle5', 'reaxff3_penalty',
123 -
    'reaxff3_angle4'
124 -
]
207 +
def read_lammps_format(lines, tolerances=None):
208 +
    """Read a reaxff file, in lammps format, to a standardised potential dictionary.
125 209
126 -
_torkeys = [
127 -
    'idx1', 'idx2', 'idx3', 'idx4', 'reaxff4_torsion1', 'reaxff4_torsion2',
128 -
    'reaxff4_torsion3', 'reaxff4_torsion4', 'reaxff4_torsion5', 'dummy1',
129 -
    'dummy2'
130 -
]
210 +
    Parameters
211 +
    ----------
212 +
    lines : list[str]
213 +
    tolerances : dict or None
214 +
        tolerances to set, that are not specified in the file.
131 215
132 -
_hbkeys = [
133 -
    'idx1', 'idx2', 'idx3', 'reaxff3_hbond1', 'reaxff3_hbond2',
134 -
    'reaxff3_hbond3', 'reaxff3_hbond4'
135 -
]
216 +
    Returns
217 +
    -------
218 +
    dict
136 219
220 +
    Notes
221 +
    -----
222 +
    Default tolerances:
137 223
138 -
def read_reaxff_file(inpath, reaxfftol=None):
139 -
    """
224 +
    - anglemin: 0.001
225 +
    - angleprod: 0.001
226 +
    - hbondmin: 0.01
227 +
    - hbonddist: 7.5
228 +
    - torsionprod: 1e-05
140 229
141 -
    :param inpath: path to reaxff file (in standard (lammps) format)
142 -
    :param reaxfftol: additional tolerance parameters
143 -
    :return:
144 230
    """
231 +
    output = {
232 +
        "description": lines[0],
233 +
        "global": {},
234 +
        "species": ["X core"],  # X is always first
235 +
        "1body": {},
236 +
        "2body": {},
237 +
        "3body": {},
238 +
        "4body": {},
239 +
    }
145 240
146 -
    toldict = _create_tol_dict(reaxfftol)
147 -
148 -
    with open(inpath, 'r') as f:
149 -
        # Descript Initial Line
150 -
        descript = f.readline()
151 -
        # Read Parameters
152 -
        # ----------
153 -
        if int(f.readline().split()[0]) != len(_paramkeys):
154 -
            raise IOError('Expecting {} general parameters'.format(
155 -
                len(_paramkeys)))
156 -
        reaxff_par = {}
157 -
        for key in _paramkeys:
158 -
            reaxff_par[key] = float(f.readline().split()[0])
159 -
160 -
        # Read Species Information
161 -
        # --------------------
162 -
        spec_dict = _read_species_info(f)
163 -
164 -
        # Read Bond Information
165 -
        # --------------------
166 -
        bond_dict = _read_bond_info(f)
167 -
168 -
        # Read Off-Diagonal Information
169 -
        # --------------------
170 -
        od_dict = _read_off_diagonal_info(f)
171 -
172 -
        # Read Angle Information
173 -
        # --------------------
174 -
        angle_dict = _read_angle_info(f)
175 -
176 -
        # Read Torsion Information
177 -
        # --------------------
178 -
        torsion_dict = _read_torsion_info(f)
179 -
180 -
        # Read HBond Information
181 -
        # --------------------
182 -
        hbond_dict = _read_hbond_info(f)
183 -
184 -
        return {
185 -
            "descript": descript.strip(),
186 -
            "tolerances": toldict,
187 -
            "params": reaxff_par,
188 -
            "species":
189 -
            spec_dict,  # spec_df.reset_index().to_dict(orient='list'),
190 -
            "bonds":
191 -
            bond_dict,  # bond_df.reset_index().to_dict(orient='list'),
192 -
            "off-diagonals":
193 -
            od_dict,  # od_df.reset_index().to_dict(orient='list'),
194 -
            "hbonds":
195 -
            hbond_dict,  # hbond_df.reset_index().to_dict(orient='list'),
196 -
            "torsions":
197 -
            torsion_dict,  # torsion_df.reset_index().to_dict(orient='list'),
198 -
            "angles":
199 -
            angle_dict,  # angle_df.reset_index().to_dict(orient='list')
200 -
        }
201 -
202 -
203 -
def _read_species_info(f):
204 -
    nspec = readnumline(f, 'species')
205 -
    skip(f, 3)
206 -
    spec_values = []
207 -
    for i in range(nspec):
208 -
        values = f.readline().split()
209 -
        symbol = values.pop(0)
210 -
        idx = 0 if symbol == 'X' else i + 1
211 -
        spec_value = [idx, symbol]
212 -
        spec_value.extend([float(v) for v in values])
213 -
        spec_value.extend(readlist(f))
214 -
        spec_value.extend(readlist(f))
215 -
        spec_value.extend(readlist(f))
216 -
        spec_values.append(spec_value)
217 -
        if len(spec_values[i]) != len(_speckeys):
241 +
    lineno = 1
242 +
243 +
    # Global parameters
244 +
    if lines[lineno].split()[0] != str(len(KEYS_GLOBAL)):
245 +
        raise IOError("Expecting {} global parameters".format(len(KEYS_GLOBAL)))
246 +
247 +
    for key in KEYS_GLOBAL:
248 +
        lineno += 1
249 +
        output["global"][key] = float(lines[lineno].split()[0])
250 +
251 +
    output["global"][
252 +
        "reaxff2_pen3"
253 +
    ] = 1.0  # this is not provided by lammps, but is used by GULP
254 +
255 +
    tolerances = tolerances or {}
256 +
    output["global"].update({k: tolerances.get(k, v) for k, v in DEFAULT_TOLERANCES})
257 +
258 +
    # one-body parameters
259 +
    lineno += 1
260 +
    num_species = int(lines[lineno].split()[0])
261 +
    lineno += 3
262 +
    idx = 1
263 +
    for i in range(num_species):
264 +
        lineno += 1
265 +
        symbol, values = lines[lineno].split(None, 1)
266 +
        if symbol == "X":
267 +
            species_idx = 0  # the X symbol is always assigned index 0
268 +
        else:
269 +
            species_idx = idx
270 +
            idx += 1
271 +
            output["species"].append(symbol + " core")
272 +
        values = split_numbers(values)
273 +
        for _ in range(3):
274 +
            lineno += 1
275 +
            values.extend(split_numbers(lines[lineno]))
276 +
277 +
        if len(values) != len(KEYS_1BODY):
218 278
            raise Exception(
219 -
                'number of values different than expected for species {0}'.
220 -
                format(symbol))
221 -
222 -
    # spec_df = pd.DataFrame(spec_values, columns=speckey).set_index('idx')
223 -
    # spec_df['reaxff1_lonepair1'] = 0.5 * (spec_df.reaxff1_valence3 - spec_df.reaxff1_valence1)
224 -
    spec_dict = {k: v for k, v in zip(_speckeys, transpose(spec_values))}
225 -
    spec_dict['reaxff1_lonepair1'] = (
226 -
        0.5 * (np.array(spec_dict["reaxff1_valence3"]) - np.array(
227 -
            spec_dict["reaxff1_valence1"]))).tolist()
228 -
    return spec_dict
229 -
230 -
231 -
def _read_bond_info(f):
232 -
    # bondcode = ['idx1', 'idx2', 'Edis1', 'Edis2', 'Edis3',
233 -
    #             'pbe1', 'pbo5', '13corr', 'pbo6', 'kov',
234 -
    #             'pbe2', 'pbo3', 'pbo4', 'nu', 'pbo1', 'pbo2',
235 -
    #             'ovcorr', 'nu']
236 -
    # bonddescript = ['idx1', 'idx2', 'Sigma-bond dissociation energy', 'Pi-bond dissociation energy',
237 -
    #                 'Double pi-bond dissociation energy',
238 -
    #                 'Bond energy', 'Double pi bond order', '1,3-Bond order correction', 'Double pi bond order',
239 -
    #                 'Overcoordination penalty',
240 -
    #                 'Bond energy', 'Pi bond order', 'Pi bond order', 'dummy', 'Sigma bond order',
241 -
    #                 'Sigma bond order',
242 -
    #                 'Overcoordination BO correction', 'dummy']
243 -
    # bond_lookup_df = pd.DataFrame(np.array([bondcode, bonddescript]).T, index=bondkey)
244 -
    nbond = readnumline(f, 'bonds')
245 -
    skip(f, 1)
246 -
    bond_values = []
247 -
    for i in range(nbond):
248 -
        values = readlist(f)
249 -
        id1 = values.pop(0)
250 -
        id2 = values.pop(0)
251 -
        bond_values.append([int(id1), int(id2)] + values + readlist(f))
252 -
        if len(bond_values[i]) != len(_bondkeys):
279 +
                "number of values different than expected for species {0}, "
280 +
                "{1} != {2}".format(symbol, len(values), len(KEYS_1BODY))
281 +
            )
282 +
283 +
        key_map = {k: v for k, v in zip(KEYS_1BODY, values)}
284 +
        key_map["reaxff1_lonepair1"] = 0.5 * (
285 +
            key_map["reaxff1_valence3"] - key_map["reaxff1_valence1"]
286 +
        )
287 +
288 +
        output["1body"][str(species_idx)] = key_map
289 +
290 +
    # two-body bond parameters
291 +
    lineno += 1
292 +
    num_lines = int(lines[lineno].split()[0])
293 +
    lineno += 2
294 +
    for _ in range(num_lines):
295 +
        values = split_numbers(lines[lineno]) + split_numbers(lines[lineno + 1])
296 +
        species_idx1 = int(values.pop(0))
297 +
        species_idx2 = int(values.pop(0))
298 +
        key_name = "{}-{}".format(species_idx1, species_idx2)
299 +
        lineno += 2
300 +
301 +
        if len(values) != len(KEYS_2BODY_BONDS):
253 302
            raise Exception(
254 -
                'number of values different than expected for bond')
255 -
    # bond_df = pd.DataFrame(bond_values, columns=bondkey).set_index(['idx1', 'idx2'])
256 -
    bond_dict = {k: v for k, v in zip(_bondkeys, transpose(bond_values))}
257 -
    return bond_dict
258 -
259 -
260 -
def _read_off_diagonal_info(f):
261 -
    nod = readnumline(f, 'off-diagonal')
262 -
    od_values = []
263 -
    for i in range(nod):
264 -
        values = readlist(f)
265 -
        id1 = int(values.pop(0))
266 -
        id2 = int(values.pop(0))
267 -
        od_values.append([id1, id2] + values)
268 -
        if len(od_values[i]) != len(_odkeys):
303 +
                "number of bond values different than expected for key {0}, "
304 +
                "{1} != {2}".format(key_name, len(values), len(KEYS_2BODY_BONDS))
305 +
            )
306 +
307 +
        output["2body"][key_name] = {k: v for k, v in zip(KEYS_2BODY_BONDS, values)}
308 +
309 +
    # two-body off-diagonal parameters
310 +
    num_lines = int(lines[lineno].split()[0])
311 +
    lineno += 1
312 +
    for _ in range(num_lines):
313 +
        values = split_numbers(lines[lineno])
314 +
        species_idx1 = int(values.pop(0))
315 +
        species_idx2 = int(values.pop(0))
316 +
        key_name = "{}-{}".format(species_idx1, species_idx2)
317 +
        lineno += 1
318 +
319 +
        if len(values) != len(KEYS_2BODY_OFFDIAG):
269 320
            raise Exception(
270 -
                'number of values different than expected for off-diagonal')
271 -
    # od_df = pd.DataFrame(od_values, columns=odkey).set_index(['idx1', 'idx2'])
272 -
    od_dict = {k: v for k, v in zip(_odkeys, transpose(od_values))}
273 -
    return od_dict
274 -
275 -
276 -
def _read_angle_info(f):
277 -
    nangle = readnumline(f, 'angle')
278 -
    angle_values = []
279 -
    for i in range(nangle):
280 -
        values = readlist(f)
281 -
        id1 = int(values.pop(0))
282 -
        id2 = int(values.pop(0))
283 -
        id3 = int(values.pop(0))
284 -
        angle_values.append([id1, id2, id3] + values)
285 -
        if len(angle_values[i]) != len(_anglekeys):
321 +
                "number of off-diagonal values different than expected for key {0} (line {1}), "
322 +
                "{2} != {3}".format(
323 +
                    key_name, lineno - 1, len(values), len(KEYS_2BODY_OFFDIAG)
324 +
                )
325 +
            )
326 +
327 +
        output["2body"].setdefault(key_name, {}).update(
328 +
            {k: v for k, v in zip(KEYS_2BODY_OFFDIAG, values)}
329 +
        )
330 +
331 +
    # three-body angle parameters
332 +
    num_lines = int(lines[lineno].split()[0])
333 +
    lineno += 1
334 +
    for _ in range(num_lines):
335 +
        values = split_numbers(lines[lineno])
336 +
        species_idx1 = int(values.pop(0))
337 +
        species_idx2 = int(values.pop(0))
338 +
        species_idx3 = int(values.pop(0))
339 +
        key_name = "{}-{}-{}".format(species_idx1, species_idx2, species_idx3)
340 +
        lineno += 1
341 +
342 +
        if len(values) != len(KEYS_3BODY_ANGLES):
286 343
            raise Exception(
287 -
                'number of values different than expected for angle')
288 -
    # angle_df = pd.DataFrame(angle_values, columns=anglekey).set_index(['idx1', 'idx2', 'idx3'])
289 -
    angle_dict = {k: v for k, v in zip(_anglekeys, transpose(angle_values))}
290 -
    return angle_dict
291 -
292 -
293 -
def _read_torsion_info(f):
294 -
    ntors = readnumline(f, 'torsion')
295 -
    torsion_values = []
296 -
    for i in range(ntors):
297 -
        values = readlist(f)
298 -
        species1 = int(values.pop(0))
299 -
        species2 = int(values.pop(0))
300 -
        species3 = int(values.pop(0))
301 -
        species4 = int(values.pop(0))
302 -
        torsion_values.append(
303 -
            [species1, species2, species3, species4] + values)
304 -
        if len(torsion_values[i]) != len(_torkeys):
344 +
                "number of angle values different than expected for key {0} (line {1}), "
345 +
                "{2} != {3}".format(
346 +
                    key_name, lineno - 1, len(values), len(KEYS_3BODY_ANGLES)
347 +
                )
348 +
            )
349 +
350 +
        output["3body"].setdefault(key_name, {}).update(
351 +
            {k: v for k, v in zip(KEYS_3BODY_ANGLES, values)}
352 +
        )
353 +
354 +
    # four-body torsion parameters
355 +
    num_lines = int(lines[lineno].split()[0])
356 +
    lineno += 1
357 +
    for _ in range(num_lines):
358 +
        values = split_numbers(lines[lineno])
359 +
        species_idx1 = int(values.pop(0))
360 +
        species_idx2 = int(values.pop(0))
361 +
        species_idx3 = int(values.pop(0))
362 +
        species_idx4 = int(values.pop(0))
363 +
        key_name = "{}-{}-{}-{}".format(
364 +
            species_idx1, species_idx2, species_idx3, species_idx4
365 +
        )
366 +
        lineno += 1
367 +
368 +
        if len(values) != len(KEYS_4BODY_TORSION):
305 369
            raise Exception(
306 -
                'number of values different than expected for torsion')
307 -
    # torsion_df = pd.DataFrame(torsion_values, columns=torkey).set_index(['idx1', 'idx2', 'idx3', 'idx4'])
308 -
    torsion_dict = {k: v for k, v in zip(_torkeys, transpose(torsion_values))}
309 -
    return torsion_dict
310 -
311 -
312 -
def _read_hbond_info(f):
313 -
    nhb = readnumline(f, 'hbond')
314 -
    hbond_values = []
315 -
    for i in range(nhb):
316 -
        values = readlist(f)
317 -
        species1 = int(values.pop(0))
318 -
        species2 = int(values.pop(0))
319 -
        species3 = int(values.pop(0))
320 -
        hbond_values.append([species1, species2, species3] + values)
321 -
        if len(hbond_values[i]) != len(_hbkeys):
370 +
                "number of torsion values different than expected for key {0} (line {1}), "
371 +
                "{2} != {3}".format(
372 +
                    key_name, lineno - 1, len(values), len(KEYS_4BODY_TORSION)
373 +
                )
374 +
            )
375 +
376 +
        output["4body"].setdefault(key_name, {}).update(
377 +
            {k: v for k, v in zip(KEYS_4BODY_TORSION, values)}
378 +
        )
379 +
380 +
    # three-body h-bond parameters
381 +
    num_lines = int(lines[lineno].split()[0])
382 +
    lineno += 1
383 +
    for _ in range(num_lines):
384 +
        values = split_numbers(lines[lineno])
385 +
        species_idx1 = int(values.pop(0))
386 +
        species_idx2 = int(values.pop(0))
387 +
        species_idx3 = int(values.pop(0))
388 +
        key_name = "{}-{}-{}".format(species_idx1, species_idx2, species_idx3)
389 +
        lineno += 1
390 +
391 +
        if len(values) != len(KEYS_3BODY_HBOND):
322 392
            raise Exception(
323 -
                'number of values different than expected for hbond {0},{1},{2}'
324 -
                .format(species1, species2, species3))
325 -
    # hbond_df = pd.DataFrame(hbond_values, columns=hbkey).set_index(['idx1', 'idx2', 'idx3'])
326 -
    hbond_dict = {k: v for k, v in zip(_hbkeys, transpose(hbond_values))}
327 -
    return hbond_dict
328 -
329 -
330 -
_tolerance_defaults = {
331 -
    "anglemin": 0.001,
332 -
    "angleprod": 0.000001,
333 -
    "hbondmin": 0.01,
334 -
    "hbonddist": 7.5,
335 -
    "torsionprod": 0.000000001
336 -
    # NB: needs to be lower to get comparable energy to lammps,
337 -
    # but then won't optimize
338 -
}
339 -
340 -
341 -
def _create_tol_dict(reaxfftol):
342 -
    reaxfftol = {} if reaxfftol is None else reaxfftol.copy()
343 -
    toldict = {}
344 -
    for key, val in _tolerance_defaults.items():
345 -
        if key in reaxfftol:
346 -
            toldict[key] = reaxfftol[key]
393 +
                "number of h-bond values different than expected for key {0} (line {1}), "
394 +
                "{2} != {3}".format(
395 +
                    key_name, lineno - 1, len(values), len(KEYS_3BODY_HBOND)
396 +
                )
397 +
            )
398 +
399 +
        output["3body"].setdefault(key_name, {}).update(
400 +
            {k: v for k, v in zip(KEYS_3BODY_HBOND, values)}
401 +
        )
402 +
403 +
    return output
404 +
405 +
406 +
def format_lammps_value(value):
407 +
    return "{:.4f}".format(value)
408 +
409 +
410 +
def write_lammps_format(data):
411 +
    """Write a reaxff file, in lammps format, from a standardised potential dictionary."""
412 +
    # validate dictionary
413 +
    validate_against_schema(data, "reaxff.schema.json")
414 +
415 +
    output = [data["description"]]
416 +
417 +
    # Global parameters
418 +
    output.append("{} ! Number of general parameters".format(len(KEYS_GLOBAL)))
419 +
    for key in KEYS_GLOBAL:
420 +
        output.append("{0:.4f} ! {1}".format(data["global"][key], key))
421 +
422 +
    # one-body parameters
423 +
    output.extend(
424 +
        [
425 +
            "{0} ! Nr of atoms; cov.r; valency;a.m;Rvdw;Evdw;gammaEEM;cov.r2;#".format(
426 +
                len(data["species"])
427 +
            ),
428 +
            "alfa;gammavdW;valency;Eunder;Eover;chiEEM;etaEEM;n.u.",
429 +
            "cov r3;Elp;Heat inc.;n.u.;n.u.;n.u.;n.u.",
430 +
            "ov/un;val1;n.u.;val3,vval4",
431 +
        ]
432 +
    )
433 +
    idx_map = {}
434 +
    i = 1
435 +
    x_species_line = None
436 +
    for idx, species in enumerate(data["species"]):
437 +
        if species.endswith("shell"):
438 +
            raise ValueError(
439 +
                "only core species can be used for reaxff, not shell: {}".format(
440 +
                    species
441 +
                )
442 +
            )
443 +
        species = species[:-5]
444 +
        # X is not always present in 1body, even if it is used in nbody terms
445 +
        # see e.g. https://github.com/lammps/lammps/blob/master/potentials/ffield.reax.cho
446 +
        if species == "X" and str(idx) not in data["1body"]:
447 +
            species_lines = []
347 448
        else:
348 -
            toldict[key] = val
349 -
    return toldict
350 -
351 -
352 -
# pylint: disable=too-many-locals
353 -
def write_lammps(data):
354 -
    """write reaxff data in LAMMPS input format
355 -
356 -
    :param data: dictionary of data
357 -
    :rtype: str
358 -
    """
359 -
    outstr = ""
360 -
361 -
    def regex(x):
362 -
        return " ".join(["{:.4f}"] * x)
363 -
364 -
    outstr += ("{}".format(data["descript"]))
365 -
    # if species_filter:
366 -
    #     outstr += ("#  (Filtered by: {})\n".format(species_filter))
367 -
    outstr += "\n"
368 -
369 -
    outstr += "{} ! Number of general parameters\n".format(len(_paramkeys))
370 -
    for key in _paramkeys:
371 -
        outstr += "{0:.4f} ! {1}\n".format(data["params"][key], key)
372 -
373 -
    outstr += '{0} ! Nr of atoms; cov.r; valency;a.m;Rvdw;Evdw;gammaEEM;cov.r2;#\n'.format(
374 -
        len(data["species"][_speckeys[0]]))
375 -
    outstr += 'alfa;gammavdW;valency;Eunder;Eover;chiEEM;etaEEM;n.u.\n'
376 -
    outstr += 'cov r3;Elp;Heat inc.;n.u.;n.u.;n.u.;n.u.\n'
377 -
    outstr += 'ov/un;val1;n.u.;val3,vval4\n'
378 -
379 -
    spec_data = transpose(
380 -
        [data['species'][key] for key in _speckeys if key != 'idx'])
381 -
382 -
    for spec in spec_data:
383 -
        outstr += spec[0] + ' ' + regex(8).format(*spec[1:9]) + '\n'
384 -
        outstr += regex(8).format(*spec[9:17]) + '\n'
385 -
        outstr += regex(8).format(*spec[17:25]) + '\n'
386 -
        outstr += regex(8).format(*spec[25:33]) + '\n'
387 -
388 -
    outstr += '{0} ! Nr of bonds; Edis1;LPpen;n.u.;pbe1;pbo5;13corr;pbo6\n'.format(
389 -
        len(data["bonds"][_bondkeys[0]]))
390 -
    outstr += 'pbe2;pbo3;pbo4;n.u.;pbo1;pbo2;ovcorr\n'
391 -
392 -
    bond_data = transpose([data['bonds'][key] for key in _bondkeys])
393 -
394 -
    for bond in bond_data:
395 -
        outstr += '{} {} '.format(
396 -
            bond[0], bond[1]) + regex(8).format(*bond[2:10]) + '\n'
397 -
        outstr += regex(8).format(*bond[10:18]) + '\n'
398 -
399 -
    outstr += '{0} ! Nr of off-diagonal terms; Ediss;Ro;gamma;rsigma;rpi;rpi2\n'.format(
400 -
        len(data["off-diagonals"][_odkeys[0]]))
401 -
402 -
    od_data = transpose([data['off-diagonals'][key] for key in _odkeys])
403 -
404 -
    for od in od_data:
405 -
        outstr += '{} {} '.format(od[0],
406 -
                                  od[1]) + regex(6).format(*od[2:8]) + '\n'
407 -
408 -
    outstr += '{0} ! Nr of angles;at1;at2;at3;Thetao,o;ka;kb;pv1;pv2\n'.format(
409 -
        len(data["angles"][_anglekeys[0]]))
410 -
411 -
    angle_data = transpose([data['angles'][key] for key in _anglekeys])
412 -
413 -
    for angle in angle_data:
414 -
        outstr += '{} {} {} '.format(*angle[0:3]) + regex(7).format(
415 -
            *angle[3:10]) + '\n'
416 -
417 -
    outstr += '{0} ! Nr of torsions;at1;at2;at3;at4;;V1;V2;V3;V2(BO);vconj;n.u;n\n'.format(
418 -
        len(data["torsions"][_torkeys[0]]))
419 -
420 -
    torsion_data = transpose([data['torsions'][key] for key in _torkeys])
421 -
422 -
    for tor in torsion_data:
423 -
        outstr += '{} {} {} {} '.format(*tor[0:4]) + regex(7).format(
424 -
            *tor[4:11]) + '\n'
425 -
426 -
    outstr += '{0} ! Nr of hydrogen bonds;at1;at2;at3;Rhb;Dehb;vhb1\n'.format(
427 -
        len(data["hbonds"][_hbkeys[0]]))
428 -
429 -
    hbond_data = transpose([data['hbonds'][key] for key in _hbkeys])
430 -
431 -
    for hbond in hbond_data:
432 -
        outstr += '{} {} {} '.format(*hbond[0:3]) + regex(4).format(
433 -
            *hbond[3:8]) + '\n'
434 -
435 -
    return outstr
436 -
437 -
438 -
class AttrDict(dict):
439 -
    def __init__(self, *args, **kwargs):
440 -
        super(AttrDict, self).__init__(*args, **kwargs)
441 -
        self.__dict__ = self
442 -
443 -
444 -
# pylint: disable=too-many-locals
445 -
def write_gulp(data, species_filter=None):
446 -
    """write reaxff data in GULP input format
447 -
448 -
    :param data: dictionary of data
449 -
    :param species_filter: list of symbols to filter
450 -
    :rtype: str
449 +
            species_lines = [
450 +
                species
451 +
                + " "
452 +
                + " ".join(
453 +
                    [
454 +
                        format_lammps_value(data["1body"][str(idx)][k])
455 +
                        for k in KEYS_1BODY[:8]
456 +
                    ]
457 +
                ),
458 +
                " ".join(
459 +
                    [
460 +
                        format_lammps_value(data["1body"][str(idx)][k])
461 +
                        for k in KEYS_1BODY[8:16]
462 +
                    ]
463 +
                ),
464 +
                " ".join(
465 +
                    [
466 +
                        format_lammps_value(data["1body"][str(idx)][k])
467 +
                        for k in KEYS_1BODY[16:24]
468 +
                    ]
469 +
                ),
470 +
                " ".join(
471 +
                    [
472 +
                        format_lammps_value(data["1body"][str(idx)][k])
473 +
                        for k in KEYS_1BODY[24:32]
474 +
                    ]
475 +
                ),
476 +
            ]
477 +
        if species == "X":
478 +
            # X is always index 0, but must be last in the species list
479 +
            idx_map[str(idx)] = "0"
480 +
            x_species_line = species_lines
481 +
        else:
482 +
            idx_map[str(idx)] = str(i)
483 +
            i += 1
484 +
            output.extend(species_lines)
485 +
    if x_species_line:
486 +
        output.extend(x_species_line)
487 +
488 +
    # two-body angle parameters
489 +
    suboutout = []
490 +
    for key in sorted(data["2body"]):
491 +
        subdata = data["2body"][key]
492 +
        if not set(subdata.keys()).issuperset(KEYS_2BODY_BONDS):
493 +
            continue
494 +
        suboutout.extend(
495 +
            [
496 +
                " ".join([idx_map[k] for k in key.split(INDEX_SEP)])
497 +
                + " "
498 +
                + " ".join(
499 +
                    [format_lammps_value(subdata[k]) for k in KEYS_2BODY_BONDS[:8]]
500 +
                ),
501 +
                " ".join(
502 +
                    [format_lammps_value(subdata[k]) for k in KEYS_2BODY_BONDS[8:16]]
503 +
                ),
504 +
            ]
505 +
        )
506 +
507 +
    output.extend(
508 +
        [
509 +
            "{0} ! Nr of bonds; Edis1;LPpen;n.u.;pbe1;pbo5;13corr;pbo6".format(
510 +
                int(len(suboutout) / 2)
511 +
            ),
512 +
            "pbe2;pbo3;pbo4;n.u.;pbo1;pbo2;ovcorr",
513 +
        ]
514 +
        + suboutout
515 +
    )
516 +
517 +
    # two-body off-diagonal parameters
518 +
    suboutout = []
519 +
    for key in sorted(data["2body"]):
520 +
        subdata = data["2body"][key]
521 +
        if not set(subdata.keys()).issuperset(KEYS_2BODY_OFFDIAG):
522 +
            continue
523 +
        suboutout.extend(
524 +
            [
525 +
                " ".join([idx_map[k] for k in key.split(INDEX_SEP)])
526 +
                + " "
527 +
                + " ".join(
528 +
                    [format_lammps_value(subdata[k]) for k in KEYS_2BODY_OFFDIAG]
529 +
                )
530 +
            ]
531 +
        )
532 +
533 +
    output.extend(
534 +
        [
535 +
            "{0} ! Nr of off-diagonal terms; Ediss;Ro;gamma;rsigma;rpi;rpi2".format(
536 +
                len(suboutout)
537 +
            )
538 +
        ]
539 +
        + suboutout
540 +
    )
541 +
542 +
    # three-body angle parameters
543 +
    suboutout = []
544 +
    for key in sorted(data["3body"]):
545 +
        subdata = data["3body"][key]
546 +
        if not set(subdata.keys()).issuperset(KEYS_3BODY_ANGLES):
547 +
            continue
548 +
        suboutout.extend(
549 +
            [
550 +
                " ".join([idx_map[k] for k in key.split(INDEX_SEP)])
551 +
                + " "
552 +
                + " ".join([format_lammps_value(subdata[k]) for k in KEYS_3BODY_ANGLES])
553 +
            ]
554 +
        )
555 +
556 +
    output.extend(
557 +
        ["{0} ! Nr of angles;at1;at2;at3;Thetao,o;ka;kb;pv1;pv2".format(len(suboutout))]
558 +
        + suboutout
559 +
    )
560 +
561 +
    # four-body torsion parameters
562 +
    suboutout = []
563 +
    for key in sorted(data["4body"]):
564 +
        subdata = data["4body"][key]
565 +
        if not set(subdata.keys()).issuperset(KEYS_4BODY_TORSION):
566 +
            continue
567 +
        suboutout.extend(
568 +
            [
569 +
                " ".join([idx_map[k] for k in key.split(INDEX_SEP)])
570 +
                + " "
571 +
                + " ".join(
572 +
                    [format_lammps_value(subdata[k]) for k in KEYS_4BODY_TORSION]
573 +
                )
574 +
            ]
575 +
        )
576 +
577 +
    output.extend(
578 +
        [
579 +
            "{0} ! Nr of torsions;at1;at2;at3;at4;;V1;V2;V3;V2(BO);vconj;n.u;n".format(
580 +
                len(suboutout)
581 +
            )
582 +
        ]
583 +
        + suboutout
584 +
    )
585 +
586 +
    # three-body h-bond parameters
587 +
    suboutout = []
588 +
    for key in sorted(data["3body"]):
589 +
        subdata = data["3body"][key]
590 +
        if not set(subdata.keys()).issuperset(KEYS_3BODY_HBOND):
591 +
            continue
592 +
        suboutout.extend(
593 +
            [
594 +
                " ".join([idx_map[k] for k in key.split(INDEX_SEP)])
595 +
                + " "
596 +
                + " ".join([format_lammps_value(subdata[k]) for k in KEYS_3BODY_HBOND])
597 +
            ]
598 +
        )
599 +
600 +
    output.extend(
601 +
        ["{0} ! Nr of hydrogen bonds;at1;at2;at3;Rhb;Dehb;vhb1".format(len(suboutout))]
602 +
        + suboutout
603 +
    )
604 +
605 +
    output.append("")
606 +
607 +
    return "\n".join(output)
608 +
609 +
610 +
def filter_by_species(data, species):
611 +
    """filter a potential dict by a subset of species
612 +
613 +
    Parameters
614 +
    ----------
615 +
    data : dict
616 +
        a potential or fitting dict
617 +
    species : list[str]
618 +
        the species to filter by
619 +
620 +
    Returns
621 +
    -------
622 +
    dict
623 +
        data filtered by species and with all species index keys re-indexed
624 +
625 +
    Raises
626 +
    ------
627 +
    KeyError
628 +
        if the data does not adhere to the potential or fitting jsonschema
629 +
    AssertionError
630 +
        if the filter set is not a subset of the available species
451 631
452 632
    """
633 +
    species = sorted(list(set(species)))
634 +
635 +
    if not set(species).issubset(data["species"]):
636 +
        raise AssertionError(
637 +
            "the filter set ({}) is not a subset of the available species ({})".format(
638 +
                set(species), set(data["species"])
639 +
            )
640 +
        )
453 641
    data = copy.deepcopy(data)
642 +
    indices = set([str(i) for i, s in enumerate(data["species"]) if s in species])
643 +
644 +
    def convert_indices(key):
645 +
        return INDEX_SEP.join(
646 +
            [str(species.index(data["species"][int(k)])) for k in key.split(INDEX_SEP)]
647 +
        )
648 +
649 +
    for key in ["1body", "2body", "3body", "4body"]:
650 +
        if key not in data:
651 +
            continue
652 +
        data[key] = {
653 +
            convert_indices(k): v
654 +
            for k, v in data[key].items()
655 +
            if indices.issuperset(k.split(INDEX_SEP))
656 +
        }
454 657
455 -
    descript = data["descript"]
456 -
    tol_par = data["tolerances"]
457 -
    reaxff_par = data["params"]
458 -
    id_sym_dict = {
459 -
        k: v
460 -
        for k, v in zip(data["species"]['idx'], data["species"]['symbol'])
461 -
    }
462 -
    spec_df = data["species"]
463 -
    spec_df['idxs'] = list(zip(spec_df['idx']))
464 -
    bond_df = data["bonds"]
465 -
    bond_df['idxs'] = list(zip(bond_df['idx1'], bond_df['idx2']))
466 -
    od_df = data["off-diagonals"]
467 -
    od_df['idxs'] = list(zip(od_df['idx1'], od_df['idx2']))
468 -
    hbond_df = data["hbonds"]
469 -
    hbond_df['idxs'] = list(zip(hbond_df['idx2'], hbond_df['idx1'],
470 -
                                hbond_df['idx3']))
471 -
    angle_df = data["angles"]
472 -
    angle_df['idxs'] = list(zip(angle_df['idx2'], angle_df['idx1'],
473 -
                                angle_df['idx3']))
474 -
    torsion_df = data["torsions"]
475 -
    torsion_df['idxs'] = list(zip(torsion_df['idx1'], torsion_df['idx2'],
476 -
                                  torsion_df['idx3'], torsion_df['idx4']))
477 -
478 -
    # If reaxff2_bo3 = 1 needs to be set to 0 for GULP since this is a dummy value
479 -
    def gulp_conv1(val):
480 -
        return 0.0 if abs(val - 1) < 1e-12 else val
481 -
482 -
    bond_df['reaxff2_bo3'] = [gulp_conv1(i)
483 -
                              for i in bond_df["reaxff2_bo3"]]
484 -
485 -
    # If reaxff2_bo(5,n) < 0 needs to be set to 0 for GULP since this is a dummy value
486 -
    def gulp_conv2(val):
487 -
        return 0.0 if val < 0.0 else val
488 -
489 -
    bond_df['reaxff2_bo5'] = [gulp_conv2(i)
490 -
                              for i in bond_df["reaxff2_bo5"]]
491 -
492 -
    # TODO, this wasn't part of the original script, and should be better understood
493 -
    # but without it, the energies greatly differ to LAMMPS (approx equal otherwise)
494 -
    def gulp_conv3(val):
495 -
        return 0.0 if val > 0.0 else val
496 -
497 -
    spec_df['reaxff1_radii3'] = [
498 -
        gulp_conv3(i) for i in spec_df['reaxff1_radii3']
499 -
    ]
500 -
501 -
    attr_dicts = _create_attr_dicts(angle_df, bond_df, hbond_df, od_df,
502 -
                                    spec_df, torsion_df)
503 -
    angle_df, bond_df, hbond_df, od_df, spec_df, torsion_df = attr_dicts  # pylint: disable=unbalanced-tuple-unpacking
504 -
505 -
    outstr = ""
506 -
507 -
    write_data = _get_write_data_func(id_sym_dict, species_filter)
508 -
509 -
    outstr = _write_header(descript, outstr, reaxff_par, species_filter, tol_par)
510 -
511 -
    outstr += "#  Species independent parameters \n"
512 -
    outstr += "#\n"
513 -
    outstr += ("reaxff0_bond     {:12.6f} {:12.6f}\n".format(
514 -
        reaxff_par['Overcoordination 1'], reaxff_par['Overcoordination 2']))
515 -
    outstr += (
516 -
        "reaxff0_over     {:12.6f} {:12.6f} {:12.6f} {:12.6f} {:12.6f}\n".
517 -
        format(reaxff_par['Valency overcoordination 2'],
518 -
               reaxff_par['Valency overcoordination 1'],
519 -
               reaxff_par['Undercoordination 1'],
520 -
               reaxff_par['Undercoordination 2'],
521 -
               reaxff_par['Undercoordination 3']))
522 -
    outstr += ("reaxff0_valence  {:12.6f} {:12.6f} {:12.6f} {:12.6f}\n".format(
523 -
        reaxff_par['Valency undercoordination'],
524 -
        reaxff_par['Valency/lone pair'], reaxff_par['Valency angle 1'],
525 -
        reaxff_par['Valency angle 2']))
526 -
    outstr += ("reaxff0_penalty  {:12.6f} {:12.6f} {:12.6f}\n".format(
527 -
        reaxff_par['Double bond/angle'],
528 -
        reaxff_par['Double bond/angle: overcoord 1'],
529 -
        reaxff_par['Double bond/angle: overcoord 2']))
530 -
    outstr += ("reaxff0_torsion  {:12.6f} {:12.6f} {:12.6f} {:12.6f}\n".format(
531 -
        reaxff_par['Torsion/BO'], reaxff_par['Torsion overcoordination 1'],
532 -
        reaxff_par['Torsion overcoordination 2'], reaxff_par['Conjugation']))
533 -
    outstr += "reaxff0_vdw      {:12.6f}\n".format(
534 -
        reaxff_par['vdWaals shielding'])
535 -
    outstr += "reaxff0_lonepair {:12.6f}\n".format(
536 -
        reaxff_par['Valency angle/lone pair'])
537 -
538 -
    outstr += "#\n"
539 -
    outstr += "#  Species parameters \n"
540 -
    outstr += "#\n"
541 -
    outstr += write_data(
542 -
        'reaxff1_radii', spec_df,
543 -
        ['reaxff1_radii1', 'reaxff1_radii2', 'reaxff1_radii3'])
544 -
545 -
    outstr += write_data(
546 -
        'reaxff1_valence', spec_df,
547 -
        ['reaxff1_valence1', 'reaxff1_valence2', 'reaxff1_valence3',
548 -
            'reaxff1_valence4'])
549 -
    outstr += write_data(
550 -
        'reaxff1_over', spec_df,
551 -
        ['reaxff1_over1', 'reaxff1_over2', 'reaxff1_over3', 'reaxff1_over4'])
552 -
    outstr += write_data('reaxff1_under kcal', spec_df, ['reaxff1_under'])
553 -
    outstr += write_data('reaxff1_lonepair kcal', spec_df,
554 -
                         ['reaxff1_lonepair1', 'reaxff1_lonepair2'])
555 -
    outstr += write_data('reaxff1_angle', spec_df,
556 -
                         ['reaxff1_angle1', 'reaxff1_angle2'])
557 -
    outstr += write_data('reaxff1_morse kcal', spec_df,
558 -
                         ['reaxff1_morse1', 'reaxff1_morse2',
559 -
                             'reaxff1_morse3', 'reaxff1_morse4'])
560 -
561 -
    outstr += "#\n"
562 -
    outstr += "# Element parameters \n"
563 -
    outstr += "#\n"
564 -
    outstr += write_data('reaxff_chi', spec_df, ['reaxff_chi'])
565 -
    outstr += write_data('reaxff_mu', spec_df, ['reaxff_mu'])
566 -
    outstr += write_data('reaxff_gamma', spec_df, ['reaxff_gamma'])
567 -
568 -
    outstr += "#\n"
569 -
    outstr += "# Bond parameters \n"
570 -
    outstr += "#\n"
571 -
    outstr += write_data(
572 -
        'reaxff2_bo over bo13', bond_df,
573 -
        ['reaxff2_bo1', 'reaxff2_bo2', 'reaxff2_bo3',
574 -
            'reaxff2_bo4', 'reaxff2_bo5', 'reaxff2_bo6'],
575 -
        conditions=lambda s: s.reaxff2_bo7 > 0.001 and s.reaxff2_bo8 > 0.001)
576 -
    outstr += write_data(
577 -
        'reaxff2_bo bo13', bond_df,
578 -
        ['reaxff2_bo1', 'reaxff2_bo2', 'reaxff2_bo3',
579 -
            'reaxff2_bo4', 'reaxff2_bo5', 'reaxff2_bo6'],
580 -
        conditions=lambda s: s.reaxff2_bo7 > 0.001 and s.reaxff2_bo8 <= 0.001)
581 -
    outstr += write_data(
582 -
        'reaxff2_bo over', bond_df,
583 -
        ['reaxff2_bo1', 'reaxff2_bo2', 'reaxff2_bo3',
584 -
            'reaxff2_bo4', 'reaxff2_bo5', 'reaxff2_bo6'],
585 -
        conditions=lambda s: s.reaxff2_bo7 <= 0.001 and s.reaxff2_bo8 > 0.001)
586 -
    outstr += write_data(
587 -
        'reaxff2_bo', bond_df,
588 -
        ['reaxff2_bo1', 'reaxff2_bo2', 'reaxff2_bo3',
589 -
            'reaxff2_bo4', 'reaxff2_bo5', 'reaxff2_bo6'],
590 -
        conditions=lambda s: s.reaxff2_bo7 <= 0.001 and s.reaxff2_bo8 <= 0.001)
591 -
    outstr += write_data(
592 -
        'reaxff2_bond kcal', bond_df,
593 -
        ['reaxff2_bond1', 'reaxff2_bond2', 'reaxff2_bond3',
594 -
            'reaxff2_bond4', 'reaxff2_bond5'])
595 -
    outstr += write_data('reaxff2_over', bond_df, ['reaxff2_over'])
596 -
    outstr += write_data(
597 -
        'reaxff2_pen kcal',
598 -
        bond_df, ['reaxff2_bo9'],
599 -
        conditions=lambda s: s.reaxff2_bo9 > 0.0,
600 -
        extra_data=[reaxff_par['Not used 1'], 1.0])
601 -
    outstr += write_data(
602 -
        'reaxff2_morse kcal', od_df,
603 -
        ['reaxff2_morse1', 'reaxff2_morse2', 'reaxff2_morse3',
604 -
            'reaxff2_morse4', 'reaxff2_morse5', 'reaxff2_morse6'])
605 -
606 -
    outstr += "#\n"
607 -
    outstr += "# Angle parameters \n"
608 -
    outstr += "#\n"
609 -
    outstr += write_data(
610 -
        'reaxff3_angle kcal', angle_df,
611 -
        ['reaxff3_angle1', 'reaxff3_angle2', 'reaxff3_angle3',
612 -
            'reaxff3_angle4', 'reaxff3_angle5'],
613 -
        conditions=lambda s: s.reaxff3_angle2 > 0.0)
614 -
    outstr += write_data('reaxff3_penalty kcal',
615 -
                         angle_df, ['reaxff3_penalty'])
616 -
    outstr += write_data(
617 -
        'reaxff3_conjugation kcal', angle_df,
618 -
        ['reaxff3_conj'],
619 -
        conditions=lambda s: abs(s.reaxff3_conj) > 1.0E-4,
620 -
        extra_data=[
621 -
            reaxff_par['Valency angle conjugation 1'],
622 -
            reaxff_par['Valency angle conjugation 3'],
623 -
            reaxff_par['Valency angle conjugation 2']
624 -
        ])
625 -
626 -
    outstr += "#\n"
627 -
    outstr += "# Hydrogen bond parameters \n"
628 -
    outstr += "#\n"
629 -
    outstr += write_data(
630 -
        'reaxff3_hbond kcal', hbond_df,
631 -
        ['reaxff3_hbond1', 'reaxff3_hbond2',
632 -
            'reaxff3_hbond3', 'reaxff3_hbond4'])
633 -
634 -
    outstr += "#\n"
635 -
    outstr += "# Torsion parameters \n"
636 -
    outstr += "#\n"
637 -
    outstr += write_data(
638 -
        'reaxff4_torsion kcal', torsion_df,
639 -
        ['reaxff4_torsion1', 'reaxff4_torsion2', 'reaxff4_torsion3',
640 -
            'reaxff4_torsion4', 'reaxff4_torsion5'])
641 -
642 -
    return outstr
643 -
644 -
645 -
def _write_header(descript, outstr, reaxff_par, species_filter, tol_par):
646 -
    outstr += "#\n"
647 -
    outstr += "#  ReaxFF force field\n"
648 -
    outstr += "#\n"
649 -
    outstr += "#  Original paper:\n"
650 -
    outstr += "#\n"
651 -
    outstr += "#  A.C.T. van Duin, S. Dasgupta, F. Lorant and W.A. Goddard III,\n"
652 -
    outstr += "#  J. Phys. Chem. A, 105, 9396-9409 (2001)\n"
653 -
    outstr += "#\n"
654 -
    outstr += "#  Parameters description:\n"
655 -
    outstr += "#\n"
656 -
    outstr += "#  {}".format(descript)
657 -
    if species_filter:
658 -
        outstr += "#  (Filtered by: {})\n".format(species_filter)
659 -
    outstr += "#\n"
660 -
    outstr += "#  Cutoffs for VDW & Coulomb terms\n"
661 -
    outstr += "#\n"
662 -
    outstr += ("reaxFFvdwcutoff {:12.4f}\n".format(
663 -
        reaxff_par['Upper Taper-radius']))
664 -
    outstr += ("reaxFFqcutoff   {:12.4f}\n".format(
665 -
        reaxff_par['Upper Taper-radius']))
666 -
    outstr += "#\n"
667 -
    outstr += "#  Bond order threshold - check anglemin as this is cutof2 given in control file\n"
668 -
    outstr += "#\n"
669 -
    outstr += (
670 -
        "reaxFFtol  {:12.10f} {:12.10f} {:12.10f} {:12.10f} {:12.10f} {:12.10f}\n"
671 -
        .format(
672 -
            0.01 * reaxff_par['bond order cutoff'], *[
673 -
                tol_par[s] for s in
674 -
                "anglemin angleprod hbondmin hbonddist torsionprod".split()
675 -
            ]))
676 -
    outstr += "#\n"
677 -
    return outstr
678 -
679 -
680 -
def _create_attr_dicts(*dicts):
681 -
    outdicts = []
682 -
    for dct in dicts:
683 -
        outdicts.append([
684 -
            AttrDict({k: v[i]
685 -
                      for k, v in dct.items()})
686 -
            for i in range(len(dct['idxs']))
687 -
        ])
688 -
    return outdicts
689 -
690 -
691 -
def _get_write_data_func(id_sym_dict, species_filter):
692 -
    def write_data(name, df, wdata, conditions=None, extra_data=None):
693 -
        if extra_data is None:
694 -
            extra_data = []
695 -
        datastr = ""
696 -
        lines = []
697 -
        for attr_dict in df:
698 -
699 -
            if conditions is not None and not conditions(attr_dict):
700 -
                continue
701 -
702 -
            try:
703 -
                symbols = [id_sym_dict[idx] for idx in attr_dict.idxs]
704 -
            except KeyError:
705 -
                raise Exception(
706 -
                    'ERROR: Species number out of bounds when getting {}'.
707 -
                    format(name))
708 -
709 -
            if species_filter and not set(symbols).issubset(species_filter):
710 -
                continue
711 -
712 -
            regex = "{:2s} core " * len(attr_dict.idxs)
713 -
            regex += "{:8.4f} " * (len(wdata) + len(extra_data)) + "\n"
714 -
            lines.append(regex.format(*symbols + [attr_dict[w]
715 -
                                                  for w in wdata] + extra_data))
716 -
717 -
        if lines:
718 -
            datastr += "{0} \n".format(name)
719 -
            for line in lines:
720 -
                datastr += line
721 -
722 -
        return datastr
723 -
724 -
    return write_data
658 +
    data["species"] = species
659 +
660 +
    return data

@@ -1,166 +1,195 @@
Loading
1 -
import numpy as np
2 1
from aiida.common.exceptions import InputValidationError
3 2
from aiida.plugins import DataFactory
3 +
import numpy as np
4 +
4 5
from aiida_lammps.calculations.lammps import BaseLammpsCalculation
5 -
from aiida_lammps.common.utils import convert_date_string, join_keywords, get_path
6 -
from aiida_lammps.validation import validate_with_json
7 -
import six
8 -
9 -
10 -
def generate_lammps_input(calc,
11 -
                          parameters,
12 -
                          potential_obj,
13 -
                          structure_filename,
14 -
                          trajectory_filename,
15 -
                          restart_filename,
16 -
                          info_filename,
17 -
                          add_thermo_keywords,
18 -
                          version_date='11 Aug 2017', **kwargs):
19 -
20 -
    pdict = parameters.get_dict()
21 -
22 -
    random_number = np.random.randint(10000000)
23 -
24 -
    # lammps_date = convert_date_string(pdict.get("lammps_version", None))
25 -
26 -
    lammps_input_file = 'units           {0}\n'.format(
27 -
        potential_obj.default_units)
28 -
    lammps_input_file += 'boundary        p p p\n'
29 -
    lammps_input_file += 'box tilt large\n'
30 -
    lammps_input_file += 'atom_style      {0}\n'.format(
31 -
        potential_obj.atom_style)
32 -
    lammps_input_file += 'read_data       {}\n'.format(structure_filename)
33 -
34 -
    lammps_input_file += potential_obj.get_input_potential_lines()
35 -
36 -
    if "neighbor" in pdict:
37 -
        lammps_input_file += "neighbor {0} {1}\n".format(pdict["neighbor"][0],
38 -
                                                         pdict["neighbor"][1])
39 -
    if "neigh_modify" in pdict:
40 -
        lammps_input_file += "neigh_modify {}\n".format(
41 -
            join_keywords(pdict["neigh_modify"]))
42 -
43 -
    # lammps_input_file += 'neighbor        0.3 bin\n'
44 -
    # lammps_input_file += 'neigh_modify    every 1 delay 0 check no\n'
45 -
46 -
    lammps_input_file += 'timestep        {}\n'.format(pdict["timestep"])
47 -
48 -
    thermo_keywords = ["step", "temp", "epair", "emol", "etotal", "press"]
49 -
    for kwd in add_thermo_keywords:
50 -
        if kwd not in thermo_keywords:
51 -
            thermo_keywords.append(kwd)
52 -
    lammps_input_file += 'thermo_style custom {}\n'.format(" ".join(thermo_keywords))
53 -
    lammps_input_file += 'thermo          1000\n'
54 -
55 -
    restart = pdict.get("restart", False)
56 -
    if restart:
57 -
        lammps_input_file += 'restart        {0} {1}\n'.format(restart, restart_filename)
58 -
59 -
    initial_temp, _, _ = get_path(pdict, ["integration", "constraints", "temp"],
60 -
                                  default=[None, None, None], raise_error=False)
61 -
62 -
    if initial_temp is not None:
63 -
        lammps_input_file += 'velocity        all create {0} {1} dist gaussian mom yes\n'.format(
64 -
            initial_temp, random_number)
65 -
        lammps_input_file += 'velocity        all scale {}\n'.format(
66 -
            initial_temp)
67 -
68 -
    lammps_input_file += 'fix             int all {0} {1} {2}\n'.format(
69 -
        get_path(pdict, ["integration", "style"]),
70 -
        join_keywords(
71 -
            get_path(pdict, ["integration", "constraints"], {}, raise_error=False)),
72 -
        join_keywords(get_path(pdict, ["integration", "keywords"], {}, raise_error=False)))
73 -
74 -
    # lammps_input_file += 'fix             int all nvt temp {0} {0} {1}\n'.format(parameters.dict.temperature,
75 -
    #                                                                              parameters.dict.thermostat_variable)
76 -
77 -
    lammps_input_file += 'run             {}\n'.format(
78 -
        parameters.dict.equilibrium_steps)
79 -
    lammps_input_file += 'reset_timestep  0\n'
80 -
81 -
    if potential_obj.atom_style == 'charge':
82 -
        dump_variables = "element x y z q"
83 -
        dump_format = "%4s  %16.10f %16.10f %16.10f %16.10f"
84 -
    else:
85 -
        dump_variables = "element x y z"
86 -
        dump_format = "%4s  %16.10f %16.10f %16.10f"
87 -
88 -
    lammps_input_file += 'dump            aiida all custom {0} {1} {2}\n'.format(
89 -
        parameters.dict.dump_rate, trajectory_filename, dump_variables)
90 -
91 -
    # TODO find exact version when changes were made
92 -
    if version_date <= convert_date_string('10 Feb 2015'):
93 -
        dump_mod_cmnd = "format"
94 -
    else:
95 -
        dump_mod_cmnd = "format line"
96 -
97 -
    lammps_input_file += 'dump_modify     aiida {0} "{1}"\n'.format(dump_mod_cmnd, dump_format)
98 -
    lammps_input_file += 'dump_modify     aiida sort id\n'
99 -
    lammps_input_file += 'dump_modify     aiida element {}\n'.format(' '.join(potential_obj.kind_elements))
100 -
101 -
    variables = pdict.get("output_variables", [])
102 -
    if variables and 'step' not in variables:
103 -
        # always include 'step', so we can sync with the `dump` data
104 -
        # NOTE `dump` includes step 0, whereas `print` starts from step 1
105 -
        variables.append('step')
106 -
    var_aliases = []
107 -
    for var in variables:
108 -
        var_alias = var.replace("[", "_").replace("]", "_")
109 -
        var_aliases.append(var_alias)
110 -
        lammps_input_file += 'variable {0} equal {1}\n'.format(var_alias, var)
111 -
    if variables:
112 -
        lammps_input_file += 'fix sys_info all print {0} "{1}" title "{2}" file {3} screen no\n'.format(
113 -
            parameters.dict.dump_rate,
114 -
            " ".join(["${{{0}}}".format(v) for v in var_aliases]),
115 -
            " ".join(var_aliases), info_filename)
116 -
117 -
    lammps_input_file += 'run             {}\n'.format(
118 -
        parameters.dict.total_steps)
119 -
120 -
    lammps_input_file += 'variable final_energy equal etotal\n'
121 -
    lammps_input_file += 'print "final_energy: ${final_energy}"\n'
122 -
123 -
    return lammps_input_file
6 +
from aiida_lammps.common.utils import convert_date_string, get_path, join_keywords
7 +
from aiida_lammps.validation import validate_against_schema
124 8
125 9
126 10
class MdCalculation(BaseLammpsCalculation):
127 -
128 -
    _generate_input_function = generate_lammps_input
129 -
130 11
    @classmethod
131 12
    def define(cls, spec):
132 13
        super(MdCalculation, cls).define(spec)
133 14
134 -
        spec.input('metadata.options.parser_name',
135 -
                   valid_type=six.string_types, default='lammps.md')
136 -
        spec.default_output_port = 'results'
137 -
138 -
        spec.output('trajectory_data',
139 -
                    valid_type=DataFactory('array.trajectory'),
140 -
                    required=True,
141 -
                    help='atomic configuration data per dump step')
142 -
        spec.output('system_data',
143 -
                    valid_type=DataFactory('array'),
144 -
                    required=False,
145 -
                    help='selected system data per dump step')
146 -
147 -
    def validate_parameters(self, param_data, potential_object):
15 +
        spec.input(
16 +
            "metadata.options.parser_name",
17 +
            valid_type=str,
18 +
            default="lammps.md",
19 +
        )
20 +
        spec.default_output_port = "results"
21 +
22 +
        spec.output(
23 +
            "trajectory_data",
24 +
            valid_type=DataFactory("lammps.trajectory"),
25 +
            required=True,
26 +
            help="atomic configuration data per dump step",
27 +
        )
28 +
        spec.output(
29 +
            "system_data",
30 +
            valid_type=DataFactory("array"),
31 +
            required=False,
32 +
            help="selected system data per dump step",
33 +
        )
34 +
35 +
    @staticmethod
36 +
    def create_main_input_content(
37 +
        parameter_data,
38 +
        potential_data,
39 +
        kind_symbols,
40 +
        structure_filename,
41 +
        trajectory_filename,
42 +
        system_filename,
43 +
        restart_filename,
44 +
    ):
45 +
46 +
        pdict = parameter_data.get_dict()
47 +
        version_date = convert_date_string(pdict.get("lammps_version", "11 Aug 2017"))
48 +
49 +
        # Geometry Setup
50 +
        lammps_input_file = "units           {0}\n".format(potential_data.default_units)
51 +
        lammps_input_file += "boundary        p p p\n"
52 +
        lammps_input_file += "box tilt large\n"
53 +
        lammps_input_file += "atom_style      {0}\n".format(potential_data.atom_style)
54 +
        lammps_input_file += "read_data       {}\n".format(structure_filename)
55 +
56 +
        # Potential Specification
57 +
        lammps_input_file += potential_data.get_input_lines(kind_symbols)
58 +
59 +
        # Pairwise neighbour list creation
60 +
        if "neighbor" in pdict:
61 +
            # neighbor skin_dist bin/nsq/multi
62 +
            lammps_input_file += "neighbor {0} {1}\n".format(
63 +
                pdict["neighbor"][0], pdict["neighbor"][1]
64 +
            )
65 +
        if "neigh_modify" in pdict:
66 +
            # e.g. 'neigh_modify every 1 delay 0 check no\n'
67 +
            lammps_input_file += "neigh_modify {}\n".format(
68 +
                join_keywords(pdict["neigh_modify"])
69 +
            )
70 +
71 +
        # Define Timestep
72 +
        lammps_input_file += "timestep        {}\n".format(pdict["timestep"])
73 +
74 +
        # Define computation/printing of thermodynamic info
75 +
        thermo_keywords = ["step", "temp", "epair", "emol", "etotal", "press"]
76 +
        for kwd in pdict.get("thermo_keywords", []):
77 +
            if kwd not in thermo_keywords:
78 +
                thermo_keywords.append(kwd)
79 +
        lammps_input_file += "thermo_style custom {}\n".format(
80 +
            " ".join(thermo_keywords)
81 +
        )
82 +
        lammps_input_file += "thermo          1000\n"
83 +
84 +
        # Define output of restart file
85 +
        restart = pdict.get("restart", False)
86 +
        if restart:
87 +
            lammps_input_file += "restart        {0} {1}\n".format(
88 +
                restart, restart_filename
89 +
            )
90 +
91 +
        # Set the initial velocities of atoms, if a temperature is set
92 +
        initial_temp, _, _ = get_path(
93 +
            pdict,
94 +
            ["integration", "constraints", "temp"],
95 +
            default=[None, None, None],
96 +
            raise_error=False,
97 +
        )
98 +
        if initial_temp is not None:
99 +
            lammps_input_file += (
100 +
                "velocity        all create {0} {1} dist gaussian mom yes\n".format(
101 +
                    initial_temp, pdict.get("rand_seed", np.random.randint(10000000))
102 +
                )
103 +
            )
104 +
            lammps_input_file += "velocity        all scale {}\n".format(initial_temp)
105 +
106 +
        # Define Equilibration Stage
107 +
        lammps_input_file += "fix             int all {0} {1} {2}\n".format(
108 +
            get_path(pdict, ["integration", "style"]),
109 +
            join_keywords(
110 +
                get_path(pdict, ["integration", "constraints"], {}, raise_error=False)
111 +
            ),
112 +
            join_keywords(
113 +
                get_path(pdict, ["integration", "keywords"], {}, raise_error=False)
114 +
            ),
115 +
        )
116 +
117 +
        lammps_input_file += "run             {}\n".format(
118 +
            parameter_data.dict.equilibrium_steps
119 +
        )
120 +
        lammps_input_file += "reset_timestep  0\n"
121 +
122 +
        if potential_data.atom_style == "charge":
123 +
            dump_variables = "element x y z q"
124 +
            dump_format = "%4s  %16.10f %16.10f %16.10f %16.10f"
125 +
        else:
126 +
            dump_variables = "element x y z"
127 +
            dump_format = "%4s  %16.10f %16.10f %16.10f"
128 +
129 +
        lammps_input_file += "dump            aiida all custom {0} {1} {2}\n".format(
130 +
            parameter_data.dict.dump_rate, trajectory_filename, dump_variables
131 +
        )
132 +
133 +
        # TODO find exact version when changes were made
134 +
        if version_date <= convert_date_string("10 Feb 2015"):
135 +
            dump_mod_cmnd = "format"
136 +
        else:
137 +
            dump_mod_cmnd = "format line"
138 +
139 +
        lammps_input_file += 'dump_modify     aiida {0} "{1}"\n'.format(
140 +
            dump_mod_cmnd, dump_format
141 +
        )
142 +
        lammps_input_file += "dump_modify     aiida sort id\n"
143 +
        lammps_input_file += "dump_modify     aiida element {}\n".format(
144 +
            " ".join(kind_symbols)
145 +
        )
146 +
147 +
        variables = pdict.get("output_variables", [])
148 +
        if variables and "step" not in variables:
149 +
            # always include 'step', so we can sync with the `dump` data
150 +
            # NOTE `dump` includes step 0, whereas `print` starts from step 1
151 +
            variables.append("step")
152 +
        var_aliases = []
153 +
        for var in variables:
154 +
            var_alias = var.replace("[", "_").replace("]", "_")
155 +
            var_aliases.append(var_alias)
156 +
            lammps_input_file += "variable {0} equal {1}\n".format(var_alias, var)
157 +
        if variables:
158 +
            lammps_input_file += 'fix sys_info all print {0} "{1}" title "{2}" file {3} screen no\n'.format(
159 +
                parameter_data.dict.dump_rate,
160 +
                " ".join(["${{{0}}}".format(v) for v in var_aliases]),
161 +
                " ".join(var_aliases),
162 +
                system_filename,
163 +
            )
164 +
165 +
        lammps_input_file += "run             {}\n".format(
166 +
            parameter_data.dict.total_steps
167 +
        )
168 +
169 +
        lammps_input_file += "variable final_energy equal etotal\n"
170 +
        lammps_input_file += 'print "final_energy: ${final_energy}"\n'
171 +
172 +
        lammps_input_file += 'print "END_OF_COMP"\n'
173 +
174 +
        return lammps_input_file
175 +
176 +
    @staticmethod
177 +
    def validate_parameters(param_data, potential_object):
148 178
        if param_data is None:
149 179
            raise InputValidationError("parameter data not set")
150 -
        validate_with_json(param_data.get_dict(), "md")
180 +
        validate_against_schema(param_data.get_dict(), "md.schema.json")
151 181
152 182
        # ensure the potential and paramters are in the same unit systems
153 183
        # TODO convert between unit systems (e.g. using https://pint.readthedocs.io)
154 -
        punits = param_data.get_dict()['units']
184 +
        punits = param_data.get_dict()["units"]
155 185
        if not punits == potential_object.default_units:
156 -
            raise InputValidationError('the units of the parameters ({}) and potential ({}) are different'.format(
157 -
                punits, potential_object.default_units
158 -
            ))
159 -
160 -
        self._retrieve_list += []
161 -
        if self.options.trajectory_name not in self._retrieve_temporary_list:
162 -
            self._retrieve_temporary_list += [self.options.trajectory_name]
163 -
        if self.options.info_filename not in self._retrieve_temporary_list:
164 -
            self._retrieve_temporary_list += [self.options.info_filename]
186 +
            raise InputValidationError(
187 +
                "the units of the parameters ({}) and potential ({}) are different".format(
188 +
                    punits, potential_object.default_units
189 +
                )
190 +
            )
165 191
166 192
        return True
193 +
194 +
    def get_retrieve_lists(self):
195 +
        return [], [self.options.trajectory_suffix, self.options.system_suffix]

@@ -1,14 +1,12 @@
Loading
1 -
from aiida.engine import CalcJob
2 -
from aiida.common.exceptions import ValidationError
3 1
from aiida.common import CalcInfo, CodeInfo
4 -
from aiida.orm import StructureData, Dict
2 +
from aiida.common.exceptions import ValidationError
3 +
from aiida.engine import CalcJob
4 +
from aiida.orm import Dict, StructureData
5 5
from aiida.plugins import DataFactory
6 -
from aiida_lammps.common.utils import convert_date_string
6 +
import numpy as np
7 +
7 8
from aiida_lammps.common.generate_structure import generate_lammps_structure
8 9
from aiida_lammps.data.potential import EmpiricalPotential
9 -
import six
10 -
11 -
import numpy as np
12 10
13 11
14 12
def get_supercell(structure, supercell_shape):
@@ -31,9 +29,9 @@
Loading
31 29
    return supercell
32 30
33 31
34 -
def get_FORCE_CONSTANTS_txt(force_constants):
32 +
def get_force_constants(force_constants):
35 33
36 -
    force_constants = force_constants.get_array('force_constants')
34 +
    force_constants = force_constants.get_array("force_constants")
37 35
38 36
    fc_shape = force_constants.shape
39 37
    fc_txt = "%4d\n" % (fc_shape[0])
@@ -49,20 +47,20 @@
Loading
49 47
def structure_to_poscar(structure):
50 48
51 49
    atom_type_unique = np.unique(
52 -
        [site.kind_name for site in structure.sites], return_index=True)[1]
50 +
        [site.kind_name for site in structure.sites], return_index=True
51 +
    )[1]
53 52
    labels = np.diff(np.append(atom_type_unique, [len(structure.sites)]))
54 53
55 -
    poscar = ' '.join(np.unique([site.kind_name for site in structure.sites]))
56 -
    poscar += '\n1.0\n'
54 +
    poscar = " ".join(np.unique([site.kind_name for site in structure.sites]))
55 +
    poscar += "\n1.0\n"
57 56
    cell = structure.cell
58 57
    for row in cell:
59 -
        poscar += '{0: 22.16f} {1: 22.16f} {2: 22.16f}\n'.format(*row)
60 -
    poscar += ' '.join(np.unique([site.kind_name for site in structure.sites])) + '\n'
61 -
    poscar += ' '.join(np.array(labels, dtype=str)) + '\n'
62 -
    poscar += 'Cartesian\n'
58 +
        poscar += "{0: 22.16f} {1: 22.16f} {2: 22.16f}\n".format(*row)
59 +
    poscar += " ".join(np.unique([site.kind_name for site in structure.sites])) + "\n"
60 +
    poscar += " ".join(np.array(labels, dtype=str)) + "\n"
61 +
    poscar += "Cartesian\n"
63 62
    for site in structure.sites:
64 -
        poscar += '{0: 22.16f} {1: 22.16f} {2: 22.16f}\n'.format(
65 -
            *site.position)
63 +
        poscar += "{0: 22.16f} {1: 22.16f} {2: 22.16f}\n".format(*site.position)
66 64
67 65
    return poscar
68 66
@@ -70,31 +68,22 @@
Loading
70 68
def parameters_to_input_file(parameters_object):
71 69
72 70
    parameters = parameters_object.get_dict()
73 -
    input_file = ('STRUCTURE FILE POSCAR\nPOSCAR\n\n')
74 -
    input_file += ('FORCE CONSTANTS\nFORCE_CONSTANTS\n\n')
75 -
    input_file += ('PRIMITIVE MATRIX\n')
76 -
    input_file += ('{} {} {} \n').format(*np.array(parameters['primitive'])[0])
77 -
    input_file += ('{} {} {} \n').format(*np.array(parameters['primitive'])[1])
78 -
    input_file += ('{} {} {} \n').format(*np.array(parameters['primitive'])[2])
79 -
    input_file += ('\n')
80 -
    input_file += ('SUPERCELL MATRIX PHONOPY\n')
81 -
    input_file += ('{} {} {} \n').format(*np.array(parameters['supercell'])[0])
82 -
    input_file += ('{} {} {} \n').format(*np.array(parameters['supercell'])[1])
83 -
    input_file += ('{} {} {} \n').format(*np.array(parameters['supercell'])[2])
84 -
    input_file += ('\n')
71 +
    input_file = "STRUCTURE FILE POSCAR\nPOSCAR\n\n"
72 +
    input_file += "FORCE CONSTANTS\nFORCE_CONSTANTS\n\n"
73 +
    input_file += "PRIMITIVE MATRIX\n"
74 +
    input_file += ("{} {} {} \n").format(*np.array(parameters["primitive"])[0])
75 +
    input_file += ("{} {} {} \n").format(*np.array(parameters["primitive"])[1])
76 +
    input_file += ("{} {} {} \n").format(*np.array(parameters["primitive"])[2])
77 +
    input_file += "\n"
78 +
    input_file += "SUPERCELL MATRIX PHONOPY\n"
79 +
    input_file += ("{} {} {} \n").format(*np.array(parameters["supercell"])[0])
80 +
    input_file += ("{} {} {} \n").format(*np.array(parameters["supercell"])[1])
81 +
    input_file += ("{} {} {} \n").format(*np.array(parameters["supercell"])[2])
82 +
    input_file += "\n"
85 83
86 84
    return input_file
87 85
88 86
89 -
def generate_LAMMPS_potential(pair_style):
90 -
91 -
    potential_file = '# Potential file generated by aiida plugin (please check citation in the orignal file)\n'
92 -
    for key, value in pair_style.dict.data.iteritems():
93 -
        potential_file += '{}    {}\n'.format(key, value)
94 -
95 -
    return potential_file
96 -
97 -
98 87
class BaseLammpsCalculation(CalcJob):
99 88
    """
100 89
    A basic plugin for calculating force constants using Lammps.
@@ -102,91 +91,147 @@
Loading
102 91
    Requirement: the node should be able to import phonopy
103 92
    """
104 93
105 -
    _INPUT_FILE_NAME = 'input.in'
106 -
    _INPUT_STRUCTURE = 'input.data'
94 +
    _INPUT_FILE_NAME = "input.in"
95 +
    _INPUT_STRUCTURE = "input.data"
107 96
108 -
    _DEFAULT_OUTPUT_FILE_NAME = 'log.lammps'
109 -
    _DEFAULT_TRAJECTORY_FILE_NAME = 'trajectory.lammpstrj'
110 -
    _DEFAULT_OUTPUT_INFO_FILE_NAME = "system_info.dump"
111 -
    _DEFAULT_OUTPUT_RESTART_FILE_NAME = 'lammps.restart'
97 +
    _DEFAULT_OUTPUT_FILE_NAME = "log.lammps"
98 +
    _DEFAULT_TRAJECTORY_FILE_NAME = "trajectory.lammpstrj"
99 +
    _DEFAULT_SYSTEM_FILE_NAME = "system_info.dump"
100 +
    _DEFAULT_RESTART_FILE_NAME = "lammps.restart"
112 101
113 -
    _retrieve_list = []
114 -
    _retrieve_temporary_list = []
115 -
    _cmdline_params = ['-in', _INPUT_FILE_NAME]
102 +
    _cmdline_params = ("-in", _INPUT_FILE_NAME)
116 103
    _stdout_name = None
117 104
118 105
    @classmethod
119 106
    def define(cls, spec):
120 107
        super(BaseLammpsCalculation, cls).define(spec)
121 -
        spec.input('structure', valid_type=StructureData, help='the structure')
122 -
        spec.input('potential', valid_type=EmpiricalPotential,
123 -
                   help='lammps potential')
124 -
        spec.input('parameters', valid_type=Dict,
125 -
                   help='the parameters', required=False)
126 -
        spec.input('metadata.options.cell_transform_filename',
127 -
                   valid_type=six.string_types, default="cell_transform.npy")
128 -
        spec.input('metadata.options.output_filename',
129 -
                   valid_type=six.string_types, default=cls._DEFAULT_OUTPUT_FILE_NAME)
130 -
        spec.input('metadata.options.trajectory_name',
131 -
                   valid_type=six.string_types, default=cls._DEFAULT_TRAJECTORY_FILE_NAME)
132 -
        spec.input('metadata.options.info_filename',
133 -
                   valid_type=six.string_types, default=cls._DEFAULT_OUTPUT_INFO_FILE_NAME)
134 -
        spec.input('metadata.options.restart_filename',
135 -
                   valid_type=six.string_types, default=cls._DEFAULT_OUTPUT_RESTART_FILE_NAME)
136 -
137 -
        spec.output('results',
138 -
                    valid_type=DataFactory('dict'),
139 -
                    required=True,
140 -
                    help='the data extracted from the main output file')
141 -
        spec.default_output_node = 'results'
142 -
143 -
        # TODO review aiidateam/aiida_core#2997, when closed, for exit code formalization
108 +
        spec.input("structure", valid_type=StructureData, help="the structure")
109 +
        spec.input("potential", valid_type=EmpiricalPotential, help="lammps potential")
110 +
        spec.input("parameters", valid_type=Dict, help="the parameters", required=False)
111 +
        spec.input(
112 +
            "metadata.options.cell_transform_filename",
113 +
            valid_type=str,
114 +
            default="cell_transform.npy",
115 +
        )
116 +
        spec.input(
117 +
            "metadata.options.output_filename",
118 +
            valid_type=str,
119 +
            default=cls._DEFAULT_OUTPUT_FILE_NAME,
120 +
        )
121 +
        spec.input(
122 +
            "metadata.options.trajectory_suffix",
123 +
            valid_type=str,
124 +
            default=cls._DEFAULT_TRAJECTORY_FILE_NAME,
125 +
        )
126 +
        spec.input(
127 +
            "metadata.options.system_suffix",
128 +
            valid_type=str,
129 +
            default=cls._DEFAULT_SYSTEM_FILE_NAME,
130 +
        )
131 +
        spec.input(
132 +
            "metadata.options.restart_filename",
133 +
            valid_type=str,
134 +
            default=cls._DEFAULT_RESTART_FILE_NAME,
135 +
        )
136 +
137 +
        spec.output(
138 +
            "results",
139 +
            valid_type=DataFactory("dict"),
140 +
            required=True,
141 +
            help="the data extracted from the main output file",
142 +
        )
143 +
        spec.default_output_node = "results"
144 144
145 145
        # Unrecoverable errors: resources like the retrieved folder or its expected contents are missing
146 146
        spec.exit_code(
147 -
            200, 'ERROR_NO_RETRIEVED_FOLDER',
148 -
            message='The retrieved folder data node could not be accessed.')
147 +
            200,
148 +
            "ERROR_NO_RETRIEVED_FOLDER",
149 +
            message="The retrieved folder data node could not be accessed.",
150 +
        )
149 151
        spec.exit_code(
150 -
            201, 'ERROR_NO_RETRIEVED_TEMP_FOLDER',
151 -
            message='The retrieved temporary folder data node could not be accessed.')
152 +
            201,
153 +
            "ERROR_NO_RETRIEVED_TEMP_FOLDER",
154 +
            message="The retrieved temporary folder data node could not be accessed.",
155 +
        )
152 156
        spec.exit_code(
153 -
            202, 'ERROR_LOG_FILE_MISSING',
154 -
            message='the main log output file was not found')
157 +
            202,
158 +
            "ERROR_LOG_FILE_MISSING",
159 +
            message="the main log output file was not found",
160 +
        )
155 161
        spec.exit_code(
156 -
            203, 'ERROR_TRAJ_FILE_MISSING',
157 -
            message='the trajectory output file was not found')
162 +
            203,
163 +
            "ERROR_TRAJ_FILE_MISSING",
164 +
            message="the trajectory output file was not found",
165 +
        )
158 166
        spec.exit_code(
159 -
            204, 'ERROR_STDOUT_FILE_MISSING',
160 -
            message='the stdout output file was not found')
167 +
            204,
168 +
            "ERROR_STDOUT_FILE_MISSING",
169 +
            message="the stdout output file was not found",
170 +
        )
161 171
        spec.exit_code(
162 -
            205, 'ERROR_STDERR_FILE_MISSING',
163 -
            message='the stderr output file was not found')
172 +
            205,
173 +
            "ERROR_STDERR_FILE_MISSING",
174 +
            message="the stderr output file was not found",
175 +
        )
164 176
165 177
        # Unrecoverable errors: required retrieved files could not be read, parsed or are otherwise incomplete
166 178
        spec.exit_code(
167 -
            300, 'ERROR_LOG_PARSING',
168 -
            message=('An error was flagged trying to parse the '
169 -
                     'main lammps output log file'))
179 +
            300,
180 +
            "ERROR_LOG_PARSING",
181 +
            message=(
182 +
                "An error was flagged trying to parse the "
183 +
                "main lammps output log file"
184 +
            ),
185 +
        )
170 186
        spec.exit_code(
171 -
            310, 'ERROR_TRAJ_PARSING',
172 -
            message=('An error was flagged trying to parse the '
173 -
                     'trajectory output file'))
187 +
            310,
188 +
            "ERROR_TRAJ_PARSING",
189 +
            message=(
190 +
                "An error was flagged trying to parse the " "trajectory output file"
191 +
            ),
192 +
        )
174 193
        spec.exit_code(
175 -
            320, 'ERROR_INFO_PARSING',
176 -
            message=('An error was flagged trying to parse the '
177 -
                     'system info output file'))
194 +
            320,
195 +
            "ERROR_INFO_PARSING",
196 +
            message=(
197 +
                "An error was flagged trying to parse the " "system info output file"
198 +
            ),
199 +
        )
178 200
179 201
        # Significant errors but calculation can be used to restart
180 202
        spec.exit_code(
181 -
            400, 'ERROR_LAMMPS_RUN',
182 -
            message='The main lammps output file flagged an error')
203 +
            400,
204 +
            "ERROR_LAMMPS_RUN",
205 +
            message="The main lammps output file flagged an error",
206 +
        )
207 +
        spec.exit_code(
208 +
            401,
209 +
            "ERROR_RUN_INCOMPLETE",
210 +
            message="The main lammps output file did not flag that the computation finished",
211 +
        )
183 212
184 -
    def validate_parameters(self, param_data, potential_object):
213 +
    @staticmethod
214 +
    def validate_parameters(param_data, potential_object):
185 215
        return True
186 216
187 217
    def prepare_extra_files(self, tempfolder, potential_object):
188 218
        return True
189 219
220 +
    def get_retrieve_lists(self):
221 +
        return [], []
222 +
223 +
    @staticmethod
224 +
    def create_main_input_content(
225 +
        parameter_data,
226 +
        potential_data,
227 +
        kind_symbols,
228 +
        structure_filename,
229 +
        trajectory_filename,
230 +
        system_filename,
231 +
        restart_filename,
232 +
    ):
233 +
        raise NotImplementedError
234 +
190 235
    def prepare_for_submission(self, tempfolder):
191 236
        """Create the input files from the input nodes passed to this instance of the `CalcJob`.
192 237
@@ -194,46 +239,49 @@
Loading
194 239
        :return: `aiida.common.CalcInfo` instance
195 240
        """
196 241
        # assert that the potential and structure have the same kind elements
197 -
        if [k.symbol for k in self.inputs.structure.kinds] != self.inputs.potential.kind_elements:
198 -
            raise ValidationError("the structure and potential are not compatible (different kind elements)")
199 -
200 -
        # Setup potential
201 -
        potential_txt = self.inputs.potential.get_potential_file()
242 +
        if self.inputs.potential.allowed_element_names is not None and not set(
243 +
            k.symbol for k in self.inputs.structure.kinds
244 +
        ).issubset(self.inputs.potential.allowed_element_names):
245 +
            raise ValidationError(
246 +
                "the structure and potential are not compatible (different kind elements)"
247 +
            )
202 248
203 249
        # Setup structure
204 250
        structure_txt, struct_transform = generate_lammps_structure(
205 -
            self.inputs.structure, self.inputs.potential.atom_style)
251 +
            self.inputs.structure, self.inputs.potential.atom_style
252 +
        )
206 253
207 -
        with open(tempfolder.get_abs_path(self.options.cell_transform_filename), 'w+b') as handle:
254 +
        with open(
255 +
            tempfolder.get_abs_path(self.options.cell_transform_filename), "w+b"
256 +
        ) as handle:
208 257
            np.save(handle, struct_transform)
209 258
210 259
        if "parameters" in self.inputs:
211 260
            parameters = self.inputs.parameters
212 261
        else:
213 262
            parameters = Dict()
214 -
        pdict = parameters.get_dict()
215 -
216 -
        # Check lammps version date in parameters
217 -
        lammps_date = convert_date_string(
218 -
            pdict.get("lammps_version", '11 Aug 2017'))
219 263
220 264
        # Setup input parameters
221 -
        input_txt = self._generate_input_function(
222 -
            parameters=parameters,
223 -
            potential_obj=self.inputs.potential,
265 +
        input_txt = self.create_main_input_content(
266 +
            parameter_data=parameters,
267 +
            potential_data=self.inputs.potential,
268 +
            kind_symbols=[kind.symbol for kind in self.inputs.structure.kinds],
224 269
            structure_filename=self._INPUT_STRUCTURE,
225 -
            trajectory_filename=self.options.trajectory_name,
226 -
            info_filename=self.options.info_filename,
270 +
            trajectory_filename=self.options.trajectory_suffix,
271 +
            system_filename=self.options.system_suffix,
227 272
            restart_filename=self.options.restart_filename,
228 -
            add_thermo_keywords=pdict.get("thermo_keywords", []),
229 -
            version_date=lammps_date)
273 +
        )
230 274
231 275
        input_filename = tempfolder.get_abs_path(self._INPUT_FILE_NAME)
232 276
233 -
        with open(input_filename, 'w') as infile:
277 +
        with open(input_filename, "w") as infile:
234 278
            infile.write(input_txt)
235 279
236 280
        self.validate_parameters(parameters, self.inputs.potential)
281 +
        retrieve_list, retrieve_temporary_list = self.get_retrieve_lists()
282 +
        retrieve_list.extend(
283 +
            [self.options.output_filename, self.options.cell_transform_filename]
284 +
        )
237 285
238 286
        # prepare extra files if needed
239 287
        self.prepare_extra_files(tempfolder, self.inputs.potential)
@@ -241,29 +289,26 @@
Loading
241 289
        # =========================== dump to file =============================
242 290
243 291
        structure_filename = tempfolder.get_abs_path(self._INPUT_STRUCTURE)
244 -
        with open(structure_filename, 'w') as infile:
292 +
        with open(structure_filename, "w") as infile:
245 293
            infile.write(structure_txt)
246 294
247 -
        if potential_txt is not None:
248 -
            potential_filename = tempfolder.get_abs_path(
249 -
                self.inputs.potential.potential_filename)
250 -
            with open(potential_filename, 'w') as infile:
251 -
                infile.write(potential_txt)
295 +
        for name, content in self.inputs.potential.get_external_files().items():
296 +
            fpath = tempfolder.get_abs_path(name)
297 +
            with open(fpath, "w") as infile:
298 +
                infile.write(content)
252 299
253 300
        # ============================ calcinfo ================================
254 301
255 302
        codeinfo = CodeInfo()
256 -
        codeinfo.cmdline_params = self._cmdline_params
303 +
        codeinfo.cmdline_params = list(self._cmdline_params)
257 304
        codeinfo.code_uuid = self.inputs.code.uuid
258 -
        codeinfo.withmpi = False  # Set lammps openmpi environment properly
305 +
        codeinfo.withmpi = self.metadata.options.withmpi
259 306
        codeinfo.stdout_name = self._stdout_name
260 307
261 308
        calcinfo = CalcInfo()
262 309
        calcinfo.uuid = self.uuid
263 -
        calcinfo.retrieve_list = self._retrieve_list + [
264 -
            self.options.output_filename,
265 -
            self.options.cell_transform_filename]
266 -
        calcinfo.retrieve_temporary_list = self._retrieve_temporary_list
310 +
        calcinfo.retrieve_list = retrieve_list
311 +
        calcinfo.retrieve_temporary_list = retrieve_temporary_list
267 312
        calcinfo.codes_info = [codeinfo]
268 313
269 314
        return calcinfo

@@ -0,0 +1,109 @@
Loading
1 +
#!/usr/bin/env python
2 +
# -*- coding: utf-8 -*-
3 +
"""Utility functions for validating JSON objects against schemas."""
4 +
import io
5 +
import json
6 +
import os
7 +
8 +
import importlib_resources
9 +
import jsonschema
10 +
11 +
from aiida_lammps.validation import schemas
12 +
13 +
14 +
def load_schema(name):
15 +
    """Read and return a JSON schema.
16 +
17 +
    If the name is an absolute path, it will be used as is, otherwise
18 +
    it will be loaded as resource from the internal json schema module.
19 +
20 +
    Parameters
21 +
    ----------
22 +
    name: str
23 +
24 +
    Returns
25 +
    -------
26 +
    dict
27 +
28 +
    """
29 +
    if os.path.isabs(name):
30 +
        with io.open(name) as jfile:
31 +
            schema = json.load(jfile)
32 +
    else:
33 +
        schema = json.loads(importlib_resources.read_text(schemas, name))
34 +
35 +
    return schema
36 +
37 +
38 +
def load_validator(schema):
39 +
    """Create a validator for a schema.
40 +
41 +
    Parameters
42 +
    ----------
43 +
    schema : str or dict
44 +
        schema or path to schema
45 +
46 +
    Returns
47 +
    -------
48 +
    jsonschema.IValidator
49 +
        the validator to use
50 +
51 +
    """
52 +
    if isinstance(schema, str):
53 +
        schema = load_schema(schema)
54 +
55 +
    validator_cls = jsonschema.validators.validator_for(schema)
56 +
    validator_cls.check_schema(schema)
57 +
58 +
    # by default, only validates lists
59 +
    def is_array(checker, instance):
60 +
        return isinstance(instance, (tuple, list))
61 +
62 +
    type_checker = validator_cls.TYPE_CHECKER.redefine("array", is_array)
63 +
    validator_cls = jsonschema.validators.extend(
64 +
        validator_cls, type_checker=type_checker
65 +
    )
66 +
67 +
    validator = validator_cls(schema=schema)
68 +
    return validator
69 +
70 +
71 +
def validate_against_schema(data, schema):
72 +
    """Validate json-type data against a schema.
73 +
74 +
    Parameters
75 +
    ----------
76 +
    data: dict
77 +
    schema: dict or str
78 +
        schema, name of schema resource, or absolute path to a schema
79 +
80 +
    Raises
81 +
    ------
82 +
    jsonschema.exceptions.SchemaError
83 +
        if the schema is invalid
84 +
    jsonschema.exceptions.ValidationError
85 +
        if the instance is invalid
86 +
87 +
    Returns
88 +
    -------
89 +
    bool
90 +
        return True if validated
91 +
92 +
93 +
    """
94 +
    validator = load_validator(schema)
95 +
    # validator.validate(data)
96 +
    errors = sorted(validator.iter_errors(data), key=lambda e: e.path)
97 +
    if errors:
98 +
        raise jsonschema.ValidationError(
99 +
            "\n".join(
100 +
                [
101 +
                    "- {} [key path: '{}']".format(
102 +
                        error.message, "/".join([str(p) for p in error.path])
103 +
                    )
104 +
                    for error in errors
105 +
                ]
106 +
            )
107 +
        )
108 +
109 +
    return True

@@ -1,81 +1,91 @@
Loading
1 1
import traceback
2 +
3 +
from aiida.orm import ArrayData, Dict
2 4
import numpy as np
3 -
from aiida.orm import Dict, TrajectoryData, ArrayData
4 5
6 +
from aiida_lammps.common.raw_parsers import convert_units, get_units_dict
7 +
from aiida_lammps.data.trajectory import LammpsTrajectory
5 8
from aiida_lammps.parsers.lammps.base import LAMMPSBaseParser
6 -
from aiida_lammps.common.raw_parsers import read_lammps_trajectory, get_units_dict
7 9
8 10
9 11
class MdParser(LAMMPSBaseParser):
10 -
    """
11 -
    Simple Parser for LAMMPS.
12 -
    """
12 +
    """Parser for LAMMPS MD calculations."""
13 13
14 14
    def __init__(self, node):
15 -
        """
16 -
        Initialize the instance of MDLammpsParser
17 -
        """
15 +
        """Initialize the instance of Lammps MD Parser."""
18 16
        super(MdParser, self).__init__(node)
19 17
20 18
    def parse(self, **kwargs):
21 -
        """
22 -
        Parses the datafolder, stores results.
23 -
        """
19 +
        """Parse the retrieved folder and store results."""
24 20
        # retrieve resources
25 -
        resources, exit_code = self.get_parsing_resources(
26 -
            kwargs, traj_in_temp=True, sys_info=True)
27 -
        if exit_code is not None:
28 -
            return exit_code
29 -
        trajectory_filename, trajectory_filepath, info_filepath = resources
21 +
        resources = self.get_parsing_resources(kwargs, traj_in_temp=True)
22 +
        if resources.exit_code is not None:
23 +
            return resources.exit_code
30 24
31 25
        # parse log file
32 26
        log_data, exit_code = self.parse_log_file()
33 27
        if exit_code is not None:
34 28
            return exit_code
35 29
36 -
        # parse trajectory file
37 -
        try:
38 -
            timestep = self.node.inputs.parameters.dict.timestep
39 -
            positions, charges, step_ids, cells, symbols, time = read_lammps_trajectory(
40 -
                trajectory_filepath, timestep=timestep,
41 -
                log_warning_func=self.logger.warning)
42 -
        except Exception:
43 -
            traceback.print_exc()
44 -
            return self.exit_codes.ERROR_TRAJ_PARSING
30 +
        traj_error = None
31 +
        if not resources.traj_paths:
32 +
            traj_error = self.exit_codes.ERROR_TRAJ_FILE_MISSING
33 +
        else:
34 +
            try:
35 +
                trajectory_data = LammpsTrajectory(resources.traj_paths[0])
36 +
                self.out("trajectory_data", trajectory_data)
37 +
            except Exception as err:
38 +
                traceback.print_exc()
39 +
                self.logger.error(str(err))
40 +
                traj_error = self.exit_codes.ERROR_TRAJ_PARSING
45 41
46 42
        # save results into node
47 43
        output_data = log_data["data"]
48 -
        if 'units_style' in output_data:
49 -
            output_data.update(get_units_dict(output_data['units_style'],
50 -
                                              ["distance", "time", "energy"]))
44 +
        if "units_style" in output_data:
45 +
            output_data.update(
46 +
                get_units_dict(
47 +
                    output_data["units_style"], ["distance", "time", "energy"]
48 +
                )
49 +
            )
51 50
        else:
52 51
            self.logger.warning("units missing in log")
53 52
        self.add_warnings_and_errors(output_data)
54 53
        self.add_standard_info(output_data)
54 +
        if "parameters" in self.node.get_incoming().all_link_labels():
55 +
            output_data["timestep_picoseconds"] = convert_units(
56 +
                self.node.inputs.parameters.dict.timestep,
57 +
                output_data["units_style"],
58 +
                "time",
59 +
                "picoseconds",
60 +
            )
55 61
        parameters_data = Dict(dict=output_data)
56 -
        self.out('results', parameters_data)
57 -
58 -
        # save trajectories into node
59 -
        trajectory_data = TrajectoryData()
60 -
        trajectory_data.set_trajectory(
61 -
            symbols, positions, stepids=step_ids, cells=cells, times=time)
62 -
        if charges is not None:
63 -
            trajectory_data.set_array('charges', charges)       
64 -
        self.out('trajectory_data', trajectory_data)
62 +
        self.out("results", parameters_data)
65 63
66 64
        # parse the system data file
67 -
        if info_filepath:
65 +
        sys_data_error = None
66 +
        if resources.sys_paths:
68 67
            sys_data = ArrayData()
69 68
            try:
70 -
                with open(info_filepath) as handle:
69 +
                with open(resources.sys_paths[0]) as handle:
71 70
                    names = handle.readline().strip().split()
72 -
                for i, col in enumerate(np.loadtxt(info_filepath, skiprows=1, unpack=True)):
71 +
                for i, col in enumerate(
72 +
                    np.loadtxt(resources.sys_paths[0], skiprows=1, unpack=True, ndmin=2)
73 +
                ):
73 74
                    sys_data.set_array(names[i], col)
74 75
            except Exception:
75 76
                traceback.print_exc()
76 -
                return self.exit_codes.ERROR_INFO_PARSING
77 -
            sys_data.set_attribute('units_style', output_data.get('units_style', None))
78 -
            self.out('system_data', sys_data)
77 +
                sys_data_error = self.exit_codes.ERROR_INFO_PARSING
78 +
            sys_data.set_attribute("units_style", output_data.get("units_style", None))
79 +
            self.out("system_data", sys_data)
79 80
80 81
        if output_data["errors"]:
81 82
            return self.exit_codes.ERROR_LAMMPS_RUN
83 +
84 +
        if traj_error:
85 +
            return traj_error
86 +
87 +
        if sys_data_error:
88 +
            return sys_data_error
89 +
90 +
        if not log_data.get("found_end", False):
91 +
            return self.exit_codes.ERROR_RUN_INCOMPLETE

@@ -1,61 +1,99 @@
Loading
1 -
from aiida.orm import Dict, ArrayData
1 +
from aiida.orm import ArrayData, Dict
2 +
import numpy as np
2 3
4 +
from aiida_lammps.common.parse_trajectory import TRAJ_BLOCK  # noqa: F401
5 +
from aiida_lammps.common.parse_trajectory import iter_trajectories
6 +
from aiida_lammps.common.raw_parsers import get_units_dict
3 7
from aiida_lammps.parsers.lammps.base import LAMMPSBaseParser
4 -
from aiida_lammps.common.raw_parsers import read_lammps_positions_and_forces_txt, get_units_dict
5 8
6 9
7 10
class ForceParser(LAMMPSBaseParser):
8 -
    """
9 -
    Simple Parser for LAMMPS.
10 -
    """
11 +
    """Parser for LAMMPS single point energy calculation."""
11 12
12 13
    def __init__(self, node):
13 -
        """
14 -
        Initialize the instance of Force LammpsParser
15 -
        """
14 +
        """Initialize the instance of Force Lammps Parser."""
16 15
        super(ForceParser, self).__init__(node)
17 16
18 17
    def parse(self, **kwargs):
19 -
        """
20 -
        Parses the datafolder, stores results.
21 -
        """
18 +
        """Parse the retrieved files and store results."""
22 19
        # retrieve resources
23 -
        resources, exit_code = self.get_parsing_resources(kwargs)
24 -
        if exit_code is not None:
25 -
            return exit_code
26 -
        trajectory_filename, trajectory_filepath, info_filepath = resources
20 +
        resources = self.get_parsing_resources(kwargs)
21 +
        if resources.exit_code is not None:
22 +
            return resources.exit_code
27 23
28 24
        # parse log file
29 25
        log_data, exit_code = self.parse_log_file()
30 26
        if exit_code is not None:
31 27
            return exit_code
32 28
33 -
        # parse trajectory file
34 -
        trajectory_txt = self.retrieved.get_object_content(trajectory_filename)
35 -
        if not trajectory_txt:
36 -
            self.logger.error("trajectory file empty")
37 -
            return self.exit_codes.ERROR_TRAJ_PARSING
38 -
        positions, forces, charges, symbols, cell2 = read_lammps_positions_and_forces_txt(
39 -
            trajectory_txt)
40 -
41 -
        # save forces and stresses into node
42 -
        array_data = ArrayData()
43 -
        array_data.set_array('forces', forces)
44 -
        if charges is not None:
45 -
            array_data.set_array('charges', charges)
46 -
        self.out('arrays', array_data)
29 +
        traj_error = None
30 +
        if not resources.traj_paths:
31 +
            traj_error = self.exit_codes.ERROR_TRAJ_FILE_MISSING
32 +
        else:
33 +
            try:
34 +
                array_data = self.parse_traj_file(resources.traj_paths[0])
35 +
                self.out("arrays", array_data)
36 +
            except IOError as err:
37 +
                self.logger.error(str(err))
38 +
                traj_error = self.exit_codes.ERROR_TRAJ_PARSING
47 39
48 40
        # save results into node
49 41
        output_data = log_data["data"]
50 -
        if 'units_style' in output_data:
51 -
            output_data.update(get_units_dict(output_data['units_style'],
52 -
                                              ["energy", "force", "distance"]))
42 +
        if "units_style" in output_data:
43 +
            output_data.update(
44 +
                get_units_dict(
45 +
                    output_data["units_style"], ["energy", "force", "distance"]
46 +
                )
47 +
            )
53 48
        else:
54 49
            self.logger.warning("units missing in log")
55 50
        self.add_warnings_and_errors(output_data)
56 51
        self.add_standard_info(output_data)
57 52
        parameters_data = Dict(dict=output_data)
58 -
        self.out('results', parameters_data)
53 +
        self.out("results", parameters_data)
59 54
60 55
        if output_data["errors"]:
61 56
            return self.exit_codes.ERROR_LAMMPS_RUN
57 +
58 +
        if traj_error:
59 +
            return traj_error
60 +
61 +
        if not log_data.get("found_end", False):
62 +
            return self.exit_codes.ERROR_RUN_INCOMPLETE
63 +
64 +
    def parse_traj_file(self, trajectory_filename):
65 +
        with self.retrieved.open(trajectory_filename, "r") as handle:
66 +
            traj_steps = list(iter_trajectories(handle))
67 +
        if not traj_steps:
68 +
            raise IOError("trajectory file empty")
69 +
        if len(traj_steps) > 1:
70 +
            raise IOError("trajectory file has multiple steps (expecting only one)")
71 +
72 +
        traj_step = traj_steps[0]  # type: TRAJ_BLOCK
73 +
74 +
        for field in ["fx", "fy", "fz"]:
75 +
            if field not in traj_step.atom_fields:
76 +
                raise IOError(
77 +
                    "trajectory file does not contain fields {}".format(field)
78 +
                )
79 +
80 +
        array_data = ArrayData()
81 +
82 +
        array_data.set_array(
83 +
            "forces",
84 +
            np.array(
85 +
                [
86 +
                    traj_step.atom_fields["fx"],
87 +
                    traj_step.atom_fields["fy"],
88 +
                    traj_step.atom_fields["fz"],
89 +
                ],
90 +
                dtype=float,
91 +
            ).T,
92 +
        )
93 +
94 +
        if "q" in traj_step.atom_fields:
95 +
            array_data.set_array(
96 +
                "charges", np.array(traj_step.atom_fields["q"], dtype=float)
97 +
            )
98 +
99 +
        return array_data

@@ -1,3 +1,4 @@
Loading
1 +
from collections.abc import Mapping
1 2
from contextlib import contextmanager
2 3
import distutils.spawn
3 4
import os
@@ -5,26 +6,32 @@
Loading
5 6
import subprocess
6 7
import sys
7 8
8 -
import six
9 -
10 9
TEST_DIR = os.path.dirname(os.path.realpath(__file__))
11 10
12 11
13 12
def lammps_version(executable="lammps"):
14 -
    """get the version of lammps
13 +
    """Get the version of lammps.
15 14
16 -
    we assume `lammps -h` returns e.g. 'LAMMPS (10 Feb 2015)' as first line
15 +
    we assume `lammps -h` returns e.g. 'LAMMPS (10 Feb 2015)' or
16 +
    'Large-scale Atomic/Molecular Massively Parallel Simulator - 5 Jun 2019'
17 17
    """
18 -
    p = subprocess.Popen([executable, "-h"], stdout=subprocess.PIPE)
19 -
    stdout, stderr = p.communicate()
20 -
    out_text = six.ensure_str(stdout)
21 -
    line = out_text.splitlines()[0]
22 -
    version = re.search(r"LAMMPS \((.*)\)", line).group(1)
23 -
    return version
18 +
    out_text = subprocess.check_output([executable, "-h"]).decode("utf8")
19 +
    match = re.search(r"LAMMPS \((.*)\)", out_text)
20 +
    if match:
21 +
        return match.group(1)
22 +
    regex = re.compile(
23 +
        r"^Large-scale Atomic/Molecular Massively Parallel Simulator - (.*)$",
24 +
        re.MULTILINE,
25 +
    )
26 +
    match = re.search(regex, out_text)
27 +
    if match:
28 +
        return match.group(1).strip()
29 +
30 +
    raise IOError("Could not find version from `{} -h`".format(executable))
24 31
25 32
26 33
def get_path_to_executable(executable):
27 -
    """ Get path to local executable.
34 +
    """Get path to local executable.
28 35
29 36
    :param executable: Name of executable in the $PATH variable
30 37
    :type executable: str
@@ -43,13 +50,12 @@
Loading
43 50
    if path is None:
44 51
        path = distutils.spawn.find_executable(executable)
45 52
    if path is None:
46 -
        raise ValueError(
47 -
            "{} executable not found in PATH.".format(executable))
53 +
        raise ValueError("{} executable not found in PATH.".format(executable))
48 54
49 55
    return os.path.abspath(path)
50 56
51 57
52 -
def get_or_create_local_computer(work_directory, name='localhost'):
58 +
def get_or_create_local_computer(work_directory, name="localhost"):
53 59
    """Retrieve or setup a local computer
54 60
55 61
    Parameters
@@ -64,20 +70,20 @@
Loading
64 70
    aiida.orm.computers.Computer
65 71
66 72
    """
67 -
    from aiida.orm import Computer
68 73
    from aiida.common import NotExistent
74 +
    from aiida.orm import Computer
69 75
70 76
    try:
71 -
        computer = Computer.objects.get(name=name)
77 +
        computer = Computer.objects.get(label=name)
72 78
    except NotExistent:
73 79
        computer = Computer(
74 -
            name=name,
75 -
            hostname='localhost',
76 -
            description=('localhost computer, '
77 -
                         'set up by aiida_lammps tests'),
78 -
            transport_type='local',
79 -
            scheduler_type='direct',
80 -
            workdir=os.path.abspath(work_directory))
80 +
            label=name,
81 +
            hostname="localhost",
82 +
            description=("localhost computer, " "set up by aiida_lammps tests"),
83 +
            transport_type="local",
84 +
            scheduler_type="direct",
85 +
            workdir=os.path.abspath(work_directory),
86 +
        )
81 87
        computer.store()
82 88
        computer.configure()
83 89
@@ -86,31 +92,34 @@
Loading
86 92
87 93
def get_or_create_code(entry_point, computer, executable, exec_path=None):
88 94
    """Setup code on localhost computer"""
89 -
    from aiida.orm import Code, Computer
90 95
    from aiida.common import NotExistent
96 +
    from aiida.orm import Code, Computer
91 97
92 -
    if isinstance(computer, six.string_types):
93 -
        computer = Computer.objects.get(name=computer)
98 +
    if isinstance(computer, str):
99 +
        computer = Computer.objects.get(label=computer)
94 100
95 101
    try:
96 102
        code = Code.objects.get(
97 -
            label='{}-{}@{}'.format(entry_point, executable, computer.name))
103 +
            label="{}-{}-{}".format(entry_point, executable, computer.label)
104 +
        )
98 105
    except NotExistent:
99 106
        if exec_path is None:
100 107
            exec_path = get_path_to_executable(executable)
101 108
        code = Code(
102 -
            input_plugin_name=entry_point,
103 -
            remote_computer_exec=[computer, exec_path],
109 +
            input_plugin_name=entry_point, remote_computer_exec=[computer, exec_path]
104 110
        )
105 -
        code.label = '{}-{}@{}'.format(
106 -
            entry_point, executable, computer.name)
111 +
        code.label = "{}-{}-{}".format(entry_point, executable, computer.label)
107 112
        code.store()
108 113
109 114
    return code
110 115
111 116
112 -
def get_default_metadata(max_num_machines=1, max_wallclock_seconds=1800, with_mpi=False,
113 -
                         num_mpiprocs_per_machine=1):
117 +
def get_default_metadata(
118 +
    max_num_machines=1,
119 +
    max_wallclock_seconds=1800,
120 +
    with_mpi=False,
121 +
    num_mpiprocs_per_machine=1,
122 +
):
114 123
    """
115 124
    Return an instance of the metadata dictionary with the minimally required parameters
116 125
    for a CalcJob and set to default values unless overridden
@@ -123,14 +132,27 @@
Loading
123 132
    :rtype: dict
124 133
    """
125 134
    return {
126 -
        'options': {
127 -
            'resources': {
128 -
                'num_machines': int(max_num_machines),
129 -
                'num_mpiprocs_per_machine': int(num_mpiprocs_per_machine)
135 +
        "options": {
136 +
            "resources": {
137 +
                "num_machines": int(max_num_machines),
138 +
                "num_mpiprocs_per_machine": int(num_mpiprocs_per_machine),
130 139
            },
131 -
            'max_wallclock_seconds': int(max_wallclock_seconds),
132 -
            'withmpi': with_mpi,
133 -
        }}
140 +
            "max_wallclock_seconds": int(max_wallclock_seconds),
141 +
            "withmpi": with_mpi,
142 +
        }
143 +
    }
144 +
145 +
146 +
def recursive_round(ob, dp, apply_lists=False):
147 +
    """ map a function on to all values of a nested dictionary """
148 +
    if isinstance(ob, Mapping):
149 +
        return {k: recursive_round(v, dp, apply_lists) for k, v in ob.items()}
150 +
    elif apply_lists and isinstance(ob, (list, tuple)):
151 +
        return [recursive_round(v, dp, apply_lists) for v in ob]
152 +
    elif isinstance(ob, float):
153 +
        return round(ob, dp)
154 +
    else:
155 +
        return ob
134 156
135 157
136 158
class AiidaTestApp(object):
@@ -161,11 +183,11 @@
Loading
161 183
        """return manager of a temporary AiiDA environment"""
162 184
        return self._environment
163 185
164 -
    def get_or_create_computer(self, name='localhost'):
186 +
    def get_or_create_computer(self, name="localhost"):
165 187
        """Setup localhost computer"""
166 188
        return get_or_create_local_computer(self.work_directory, name)
167 189
168 -
    def get_or_create_code(self, entry_point, computer_name='localhost'):
190 +
    def get_or_create_code(self, entry_point, computer_name="localhost"):
169 191
        """Setup code on localhost computer"""
170 192
171 193
        computer = self.get_or_create_computer(computer_name)
@@ -175,12 +197,16 @@
Loading
175 197
        except KeyError:
176 198
            raise KeyError(
177 199
                "Entry point {} not recognized. Allowed values: {}".format(
178 -
                    entry_point, self._executables.keys()))
200 +
                    entry_point, self._executables.keys()
201 +
                )
202 +
            )
179 203
180 204
        return get_or_create_code(entry_point, computer, executable)
181 205
182 206
    @staticmethod
183 -
    def get_default_metadata(max_num_machines=1, max_wallclock_seconds=1800, with_mpi=False):
207 +
    def get_default_metadata(
208 +
        max_num_machines=1, max_wallclock_seconds=1800, with_mpi=False
209 +
    ):
184 210
        return get_default_metadata(max_num_machines, max_wallclock_seconds, with_mpi)
185 211
186 212
    @staticmethod
@@ -198,6 +224,7 @@
Loading
198 224
199 225
        """
200 226
        from aiida.plugins import ParserFactory
227 +
201 228
        return ParserFactory(entry_point_name)
202 229
203 230
    @staticmethod
@@ -215,6 +242,7 @@
Loading
215 242
216 243
        """
217 244
        from aiida.plugins import DataFactory
245 +
218 246
        return DataFactory(entry_point_name)(**kwargs)
219 247
220 248
    @staticmethod
@@ -228,10 +256,12 @@
Loading
228 256
229 257
        """
230 258
        from aiida.plugins import CalculationFactory
259 +
231 260
        return CalculationFactory(entry_point_name)
232 261
233 -
    def generate_calcjob_node(self, entry_point_name, retrieved,
234 -
                              computer_name='localhost', attributes=None):
262 +
    def generate_calcjob_node(
263 +
        self, entry_point_name, retrieved, computer_name="localhost", attributes=None
264 +
    ):
235 265
        """Fixture to generate a mock `CalcJobNode` for testing parsers.
236 266
237 267
        Parameters
@@ -257,25 +287,25 @@
Loading
257 287
258 288
        process = self.get_calc_cls(entry_point_name)
259 289
        computer = self.get_or_create_computer(computer_name)
260 -
        entry_point = format_entry_point_string(
261 -
            'aiida.calculations', entry_point_name)
290 +
        entry_point = format_entry_point_string("aiida.calculations", entry_point_name)
262 291
263 292
        node = CalcJobNode(computer=computer, process_type=entry_point)
264 -
        spec_options = process.spec().inputs['metadata']['options']
265 -
        # TODO post v1.0.0b2, this can be replaced with process.spec_options
266 -
        node.set_options({
267 -
            k: v.default