Lammps not sorting atoms correctly

[Leimeroth]

I am trying to run some reaxff calculations with a potential DF made by this function

def REAXFF_DF(
    path,
    pot_file_name,
    elements,
    keywords='NULL checkqeq yes',
):
    elements_str = ""
    for e in elements:
        elements_str += f" {e}"
    
    if pot_file_name is None:
        fname = path
        pot_file_name = path.split("/")[-1]
    else:
        fname = f"{path}/{pot_file_name}"

    config = [
        f"pair_style reaxff {keywords}\n",
        f"pair_coeff * * {pot_file_name} {elements_str}\n",
        f'fix qeqreaxff all qeq/reaxff 1 0.0 10.0 1e-6 reaxff\n',
    ]
    files = [fname]
    
    pot = pot = pd.DataFrame(
        {
            "Name": "REAXFF",
            "Filename": [files],
            "Model": ["Custom"],
            "Species": [elements],
            "Config": [config],
        }
    )

    return pot

This requires to set

job.input.control['atom_style'] = 'charge'

Here I noticed that lammps does not follow the order of elements given by the potential when writing the structure.inp file.
This seems like a very critical bug and used to work as far as I can remember.
In the structure_full function there is a lot of code doing the sorting:

sorted_species_list = np.array(self._potential.get_element_lst())

for species in self.structure.species:
            species_name = species.Abbreviation
            q_dict[species_name] = self.potential.get_charge(species_name)
            species_translate_list.append(self.potential.get_element_id(species_name))
        sorted_species_list = np.array(self._potential.get_element_lst())

This is completely lacking in structure_atomic for example. Instead there is
something like

        for id_eam, el_eam in enumerate(self._el_eam_lst):
            if el_eam in el_struct_lst:
                id_el = list(el_struct_lst).index(el_eam)
                el = el_obj_lst[id_el]
                el_dict[el] = id_eam + 1

which however, seems to not work.

Tagging you @pmrv and @jan-janssen, maybe one of you knows whats going on?

[jan-janssen]

@Leimeroth For traditional potentials (EAM, MEAM) it is working fine for me, can you confirm that those are also working for you? The important part is that the elements in Species are sorted in the way how they are implemented in the potential. Here is an example for SiO:

pot = pd.DataFrame({
    "Filename": [['potential_LAMMPS/1994--Nakano-A--Si-O--LAMMPS--ipr1/SiO.1994.vashishta']],
    "Model": ["NISTiprpy"],
    "Name": ["1994--Nakano-A--Si-O--LAMMPS--ipr1"],
    "Species": [['Si', 'O']],
    "Config": [['pair_style vashishta\n', 'pair_coeff * * SiO.1994.vashishta Si O\n']],
})

[Leimeroth]

So I did some tests and apparently the error occurs only when setting

job.input.control['atom_style'] = 'charge'

Seems like a pretty evil trap to me that can quickly mess up your calculations. Is there a reason why the sorting is different between different atom_styles anyway?
If not, I guess it would be good to refactor the structure writing to have the same function doing the sorting based on the Species field of the potential dataframe all the time.

[Leimeroth]

short reminder :)

[jan-janssen]

Hi @Leimeroth , given the current conferences and that none of the current pyiron developers is using the charge atom style in LAMMPS I can not promise a quick solution. So anything helps. If you can provide a simple unit test which currently does not work that would already be helpful and makes it easier for us as pyiron developers to take a look on how to fix this issue.

[jan-janssen]

@samwaseda Do you have time to take a look at this?