readcon-core 0.22.0

An oxidized single and multiple CON file reader and writer with FFI bindings for ergonomic C/C++ usage. Implements CON spec v3.
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266

"""ASE readcon-con and readcon-convel formats."""

import sys
from pathlib import Path

import numpy as np
import pytest
import readcon
from ase import Atoms
from ase.constraints import FixAtoms, FixCartesian
from ase.io import iread, read, write

_PLUGIN = Path(__file__).resolve().parents[2] / "ase-readcon"
if str(_PLUGIN) not in sys.path:
    sys.path.insert(0, str(_PLUGIN))


def _ensure_formats() -> None:
    """Register the plugin when ASE has not loaded the entry point."""
    from ase.io.formats import define_io_format, ioformats
    from ase_readcon.io import readcon_con_format, readcon_convel_format

    for name, fmt in (
        ("readcon-con", readcon_con_format),
        ("readcon-convel", readcon_convel_format),
    ):
        if name in ioformats:
            continue
        define_io_format(name, **fmt._asdict(), external=True)


_ensure_formats()


def _atoms(shift=0.0, velocity=False, cartesian=False):
    atoms = Atoms(
        "CuH",
        positions=[[0.2 + shift, 0.3, 0.4], [1.5, 1.6, 1.7]],
        cell=[4.0, 5.0, 6.0],
        pbc=[True, True, False],
    )
    atoms.set_masses([63.546, 1.008])
    if cartesian:
        atoms.set_constraint(FixCartesian(1, mask=[True, False, True]))
    else:
        atoms.set_constraint(FixAtoms(indices=[0]))
    atoms.info["label"] = "grid"
    if velocity:
        atoms.set_velocities([[0.1, -0.2, 0.3], [0.0, 0.4, -0.5]])
    return atoms


def _fixed_indices(atoms):
    found = []
    for constraint in atoms.constraints:
        if isinstance(constraint, FixAtoms):
            found.extend(int(i) for i in constraint.index)
    return found


def _cartesian(atoms):
    found = []
    for constraint in atoms.constraints:
        if isinstance(constraint, FixCartesian):
            found.append(
                (
                    int(np.asarray(constraint.index).reshape(-1)[0]),
                    [bool(v) for v in np.asarray(constraint.mask).reshape(-1)],
                )
            )
    return found


def _same_atoms(got, expect):
    assert got.get_chemical_symbols() == expect.get_chemical_symbols()
    np.testing.assert_allclose(got.positions, expect.positions)
    np.testing.assert_allclose(got.get_masses(), expect.get_masses())
    np.testing.assert_allclose(got.cell.array, expect.cell.array)
    assert list(got.pbc) == list(expect.pbc)
    assert _fixed_indices(got) == _fixed_indices(expect)
    assert _cartesian(got) == _cartesian(expect)
    assert dict(got.info) == dict(expect.info)
    np.testing.assert_allclose(got.get_velocities(), expect.get_velocities())


def test_matches_eon_reader_on_fixture(tmp_path):
    images = [_atoms(0.0), _atoms(0.25)]
    path = tmp_path / "fixture.con"
    write(path, images, format="eon")
    via_eon = read(path, index=":", format="eon")
    via_plugin = read(path, index=":", format="readcon-con")
    assert len(via_plugin) == len(via_eon)
    for got, expect in zip(via_plugin, via_eon):
        _same_atoms(got, expect)


def test_round_trip_against_ase_eon_writer(tmp_path):
    images = [_atoms(0.0), _atoms(0.25)]
    path = tmp_path / "ase.con"
    write(path, images, format="eon")
    via_eon = read(path, index=":", format="eon")
    via_plugin = read(path, index=":", format="readcon-con")
    assert len(via_plugin) == len(via_eon) == 2
    for got, expect in zip(via_plugin, via_eon):
        assert got.get_chemical_symbols() == expect.get_chemical_symbols()
        np.testing.assert_allclose(got.positions, expect.positions)
        np.testing.assert_allclose(got.get_masses(), expect.get_masses())
        np.testing.assert_allclose(got.cell.lengths(), expect.cell.lengths())
        assert _fixed_indices(got) == _fixed_indices(expect)

    out = tmp_path / "plugin.con"
    write(out, via_plugin, format="readcon-con")
    back = read(out, index=":", format="eon")
    for got, expect in zip(back, via_eon):
        np.testing.assert_allclose(got.positions, expect.positions, atol=1e-8)
        assert _fixed_indices(got) == _fixed_indices(expect)


def test_convel_velocities_and_metadata(tmp_path):
    atoms = _atoms(velocity=True)
    path = tmp_path / "v.convel"
    write(path, atoms, format="readcon-convel")
    got = read(path, format="readcon-convel")
    np.testing.assert_allclose(got.get_velocities(), atoms.get_velocities())
    assert got.info["label"] == "grid"
    assert list(got.pbc) == [True, True, False]
    con_path = tmp_path / "v.con"
    write(con_path, atoms, format="readcon-con")
    bare = read(con_path, format="readcon-con")
    velocities = bare.get_velocities()
    assert velocities is None or np.max(np.abs(velocities)) == 0.0


def test_constraints_round_trip(tmp_path):
    atoms = _atoms(cartesian=True)
    path = tmp_path / "fix.con"
    write(path, atoms, format="readcon-con")
    got = read(path, format="readcon-con")
    assert _cartesian(got) == [(1, [True, False, True])]
    np.testing.assert_allclose(got.positions, atoms.positions)
    np.testing.assert_allclose(got.get_masses(), atoms.get_masses())


def test_atom_id_round_trip_on_the_bulk_path(tmp_path):
    """Column 5 survives a single-species read of every frame and the write back."""
    images = []
    for shift, ident in ((0.0, [4, 9, 1]), (0.2, [8, 3, 6])):
        atoms = Atoms(
            "Cu3",
            positions=[[0.0 + shift, 0.0, 0.0], [1.0, 0.0, 0.0], [2.0, 0.0, 0.0]],
            cell=[5.0, 5.0, 5.0],
            pbc=True,
        )
        atoms.set_array("atom_id", np.array(ident, dtype=np.uint64))
        images.append(atoms)
    path = tmp_path / "ids.con"
    write(path, images, format="readcon-con")
    got = read(path, index=":", format="readcon-con")
    assert len(got) == 2
    np.testing.assert_array_equal(got[0].arrays["atom_id"], [4, 9, 1])
    np.testing.assert_array_equal(got[1].arrays["atom_id"], [8, 3, 6])
    np.testing.assert_allclose(got[1].positions, images[1].positions)


def test_claisen_species_order_round_trip(tmp_path):
    """Oxygen-first product keeps its masses when the reactant was carbon-first.

    The bulk reader used to stamp frame 0's elements onto later frames, so a
    carbon was built with the oxygen mass and the write raised MassMismatch.
    """
    _tests = Path(__file__).resolve().parent
    if str(_tests) not in sys.path:
        sys.path.insert(0, str(_tests))
    from test_atom_correspondence import _claisen, _frame

    reactant_rows, product_rows, *_rest = _claisen()
    reactant = _frame(reactant_rows)
    product = _frame(product_rows)
    source = tmp_path / "band.con"
    readcon.write_con(str(source), [reactant, product])
    images = read(source, index=":", format="readcon-con")
    out = tmp_path / "back.con"
    write(out, images, format="readcon-con")
    back = read(out, index=":", format="readcon-con")
    original = readcon.read_con(str(source))
    assert len(back) == len(original) == 2
    for got, frame in zip(back, original):
        assert got.get_chemical_symbols() == [atom.symbol for atom in frame.atoms]
        np.testing.assert_allclose(
            got.get_masses(), [atom.mass for atom in frame.atoms]
        )
        np.testing.assert_array_equal(
            got.arrays["atom_id"], [atom.atom_id for atom in frame.atoms]
        )
        np.testing.assert_allclose(
            got.positions, [(atom.x, atom.y, atom.z) for atom in frame.atoms]
        )


def test_oxygen_first_velocity_reaches_ase(tmp_path):
    """A velocity on a later atom is not dropped because the first atom is still."""
    _tests = Path(__file__).resolve().parent
    if str(_tests) not in sys.path:
        sys.path.insert(0, str(_tests))
    from test_atom_correspondence import _claisen, _frame

    _reactant, product_rows, *_rest = _claisen()
    product = _frame(product_rows, velocity=(5, (0.2, -0.1, 0.05)))
    atoms = product.to_ase()
    carbon = next(i for i, atom in enumerate(product.atoms) if atom.atom_id == 5)
    np.testing.assert_allclose(atoms.get_velocities()[carbon], [0.2, -0.1, 0.05])
    path = tmp_path / "product.convel"
    write(path, atoms, format="readcon-convel")
    got = read(path, format="readcon-convel")
    np.testing.assert_allclose(got.get_velocities()[carbon], [0.2, -0.1, 0.05])
    assert got.get_chemical_symbols()[0] == "O"
    np.testing.assert_allclose(got.get_masses()[0], 15.999)


def test_pt100_adatom_fixed_in_z_only(tmp_path):
    """Per-axis FixCartesian on the Pt(100) adatom survives the plugin."""
    from ase.build import add_adsorbate, fcc100

    slab = fcc100("Pt", (3, 3, 3), vacuum=7.0)
    a0 = slab.cell[0, 0] / 3
    add_adsorbate(slab, "Pt", 1.5, position=(0.5 * a0, 0.5 * a0))
    n = len(slab)
    bottom = [i for i, tag in enumerate(slab.get_tags()) if tag == 3]
    slab.set_constraint(
        [
            FixAtoms(indices=bottom),
            FixCartesian(n - 1, mask=(False, False, True)),
        ]
    )
    path = tmp_path / "slab.con"
    write(path, slab, format="readcon-con")
    got = read(path, index=":", format="readcon-con")
    assert len(got) == 1
    assert _fixed_indices(got[0]) == bottom
    assert _cartesian(got[0]) == [(n - 1, [False, False, True])]
    np.testing.assert_allclose(got[0].positions, slab.positions)


def test_multiframe_index_does_not_call_read_con(tmp_path, monkeypatch):
    images = [_atoms(i * 0.1) for i in range(5)]
    path = tmp_path / "many.con"
    write(path, images, format="readcon-con")

    def fail_read_con(*_args, **_kwargs):
        raise AssertionError("read_con loaded the whole file")

    monkeypatch.setattr(readcon, "read_con", fail_read_con)
    first = read(path, index=0, format="readcon-con")
    second = next(iread(path, index=1, format="readcon-con"))
    last = read(path, index=-1, format="readcon-con")
    middle = read(path, index="1:4:2", format="readcon-con")
    np.testing.assert_allclose(first.positions, images[0].positions)
    np.testing.assert_allclose(second.positions, images[1].positions)
    np.testing.assert_allclose(last.positions, images[-1].positions)
    assert len(middle) == 2
    np.testing.assert_allclose(middle[0].positions, images[1].positions)
    np.testing.assert_allclose(middle[1].positions, images[3].positions)
    streamed = list(iread(path, format="readcon-con"))
    assert len(streamed) == 5
    with pytest.raises(AssertionError, match="whole file"):
        readcon.read_con(str(path))