Setting selective dynamics is confusing

[Leimeroth]

Summary

I tried to work with Lammps and selective dynamics and this basically confused me how that tag is supposed to be set, because the way that seemed the most logic to me did not work. My goal was to keep a sphere inside a bulk fixed. Using a boolean array shaped (n,3) "tags" and setting s.add_tag(selective_dynamics=tags) either resulted in all atoms moving or on some weird occasions (maybe when I tried to explicitly call _set_selective_dynamics?) in no atoms moving. What I than found to be working was to create a boolean array shaped (n) and setting s.selective_dynamics[tags] = [False,False,False]. The first one however seems a lot more intuitive to me and also allows for complete flexibility which tags are set for which atom, while it seems difficult/impossible for the 2nd version.

MWE

from pyiron import Project
import numpy as np

pr = Project("TestSelDyn")

def get_center(s):
    return np.average(s.positions, axis=0)

def dist_from_center(s):
    c = get_center(s)
    dists = np.linalg.norm(s.positions-c, axis=1) 
    return dists

def radius_filter(s, r):
    dists = dist_from_center(s)
    return dists < r

def radius_selective_dyn(s, r):
    rf = radius_filter(s,r)
    s.add_tag(selective_dynamics=[True,True,True])
    s.selective_dynamics[rf] = [False,False,False]
    return rf

def radius_selective_dyn(s, r):
    rf = radius_filter(s,r)
    tags = np.array([rf, rf, rf]).T
    s.add_tag(selective_dynamics=tags)
    return tags

s = pr.create.structure.ase.bulk("Cu", cubic=True).repeat(3)
lmp = pr.create.job.Lammps("Lammps", delete_existing_job=True)
tags = radius_selective_dyn(s, 5)
lmp.structure = s
#lmp._set_selective_dynamics()
lmp.calc_md(n_ionic_steps=5000, temperature=300, pressure=0)
lmp.run()

Copy and paste for a simulation where all atoms move, despite setting selective_dynamics, switch around order of radius_selective_dyn functions to obtain one where central atoms are locked.

[jan-janssen]

As I recently had a lengthy discussion about selective dynamics with @ligerzero-ai I just want to link this here https://github.com/pyiron/pyiron_atomistics/issues/932

[ligerzero-ai]

Hi Niklas,

I have to admit that I do not use ASE for the bulk of my interactions with structure manipulation etc. specifically for the reasons that you have listed here (and other bombs lurking in their code/general user friendliness), but also some historical ones (i.e. I've got much more accumulated experience with the pymatgen manipulations compared to Atoms). But if you are willing to consider learning some additional structural manipulation syntax, and a couple import statements, I suggest simply converting it into a pymatgen structure via:

from pyiron_atomistics.atomistics.structure.atoms import pyiron_to_pymatgen, pymatgen_to_pyiron
from pymatgen.core import Structure

and then setting:

struct.add_site_property("selective_dynamics", values=[[False, False, False], [False, False, True], [False, False, False], [True, True, False]])

struct_calc = pymatgen_to_pyiron(struct)
job.structure = struct_calc

The syntax that you are suggesting is something that pymatgen supports and I like it as the more pythonic and intuitive solution. Having said this, I am reluctant to spend dev time on something that I know already has a relatively simple workaround. Of course, if there is enough demand, I guess you could pester the ASE developers to fix the clunky interactions with setting interactions.

[pmrv]

I'm afraid the clunkiness of the selective_dynamics on us and the weird SparseArray, SparseList constructions that are used for it in the back and are mostly unrelated to ASE. For this specific issue would be relatively straightforward in SparseList.add_tag.

[pmrv]

Even after #981 is merged, we should think a bit on how to make that feature more discoverable. Maybe a method that adds the selective dynamics tag and explains the semantics?

[jan-janssen]

To be provocative - do we need a separate implementation from ASE https://wiki.fysik.dtu.dk/ase/ase/constraints.html

[pmrv]

The implementation doesn't matter to me, but I suppose we want some thin interface so people do not need to import from all over them place and there's some level of discoverability.