Selective dynamics for LAMMPS seem to be broken

[jan-janssen]

I presume this this related to #981

Both setting the constraints using ASE:

_, _, z = structure.scaled_pos_xyz()
c = FixAtoms(indices=[atom.index for selector, atom in zip(z < 0.5, structure) if selector])
structure.set_constraint(c)

is broken as well as setting the constraints using pyiron:

structure.add_tag(selective_dynamics=None)
_, _, z = structure.scaled_pos_xyz()
for selector, ind in zip(z < 0.5, range(len(structure))):
    if selector:
        structure.selective_dynamics[ind] = np.array([True, True, True])
    else:
        structure.selective_dynamics[ind] = np.array([False, False, False])

Full example:

from ase.constraints import FixAtoms
from pyiron_atomistics import Project
pr = Project("test")
structure = pr.create.structure.ase.bulk("Al", cubic=True).repeat([2,2,4])
_, _, z = structure.scaled_pos_xyz()
c = FixAtoms(indices=[atom.index for selector, atom in zip(z < 0.5, structure) if selector])
structure.set_constraint(c)
job = pr.create.job.Lammps("lmp")
job.structure = structure
job.calc_md(temperature=500.0, n_ionic_steps=10000)
job.run()

The error happens in:

File ~/PycharmProjects/pyiron_atomistics/pyiron_atomistics/lammps/base.py:988, in LammpsBase._set_selective_dynamics(self)
    986 sel_dyn = np.logical_not(self.structure.selective_dynamics)
    987 # Enter loop only if constraints present
--> 988 if len(np.argwhere(np.any(sel_dyn, axis=1)).flatten()) != 0:
    989     all_indices = np.arange(len(self.structure), dtype=int)
    990     constraint_xyz = np.argwhere(np.all(sel_dyn, axis=1)).flatten()

But I presume the same happens in VASP as well. So my suggestion would be to revert #981 until we have a fix for the issue.

[pmrv]

I'm somewhat sure that it still works when setting the selective_dynamics tag on the Atoms object, which I thought was the "pyiron" way. Compatibility with the ASE constraints would be great though.

[jan-janssen]

@pmrv for me it does not work, I tired:

structure.add_tag(selective_dynamics=None)
_, _, z = structure.scaled_pos_xyz()
for selector, ind in zip(z < 0.5, range(len(structure))):
    if selector:
        structure.selective_dynamics[ind] = np.array([True, True, True])
    else:
        structure.selective_dynamics[ind] = np.array([False, False, False])

This is the code we use in the melting point script https://github.com/pyiron/pyiron_atomistics/blob/main/pyiron_atomistics/thermodynamics/interfacemethod.py#L31