Lammps interactive bugged for some structures

[Leimeroth]

Summary

Lammps modal and interactive mode calculate completely different energies for some structures. I have absolutely no idea what is happening. Example structure is attached in vasp POSCAR format. Tested for different potentials.

pyiron Version and Platform

latest git versions, Linux

Expected Behavior

Same energies, at least within numerical accuracy

Actual Behavior

0.6 eV/atom difference

Steps to Reproduce MWE with attached POSCAR:

s = pr.create.structure.read("POSCAR")

lmp = pr.create.job.Lammps("Relax", delete_existing_job=True) lmp.structure = s #lmp.potential = pot lmp.calc_minimize(pressure=0) lmp.run() print(lmp["output/generic/energy_pot"])

lmp2 = pr.create.job.Lammps("RelaxInteractive", delete_existing_job=True) lmp2.server.run_mode.interactive = True lmp2.structure = s #lmp.potential = pot lmp2.calc_minimize(pressure=0) lmp2.run() lmp2.interactive_energy_pot_getter()

POSCAR.txt

[jan-janssen]

I was able to track it down to https://github.com/pyiron/pyiron_atomistics/commit/4cd079d8b8e8ffd0032ae12afd6af484d1c2efba @liamhuber Can you take a look?

The difference is basically:

# interactive job
UnfoldingPrism(structure.cell)

# basic job
UnfoldingPrism(UnfoldingPrism(structure.cell).A)

All this happens in structure_to_lammps().

[jan-janssen]

The quick solution is to always apply:

def structure_to_lammps(structure):
    """Converts a structure to the Lammps coordinate frame"""
    prism = UnfoldingPrism(structure.cell)
    lammps_structure = structure.copy()
    lammps_structure.cell = prism.A
    lammps_structure.positions = np.matmul(structure.positions, prism.R)
    return lammps_structure

Before setting the structure.

[samwaseda]

What do you think of using something like structure._voigt for all things in LAMMPS consistently? I propose the following class:

import numpy as np
from pyiron_atomistics.atomistics.structure.atoms import Atoms as PyironAtoms
from pyiron_atomistics.lammps.structure import UnfoldingPrism


class Atoms(PyironAtoms):
    def __init__(self, *args, **kwargs):
        super().__init__(*args, **kwargs)
        self._voigt = Voigt(self)


class Voigt:
    def __init__(self, ref_structure):
        self._ref_structure = ref_structure
        self._ref_cell = ref_structure.cell.copy()
        self._prism = None

    @property
    def prism(self):
        if self._prism is None or not np.isclose(
            self._ref_cell, self.ref_structure.cell
        ).all():
            self._prism = UnfoldingPrism(self.ref_structure.cell, digits=15)
            self._ref_cell = self.ref_structure.cell.copy()
        return self._prism

    @property
    def ref_structure(self):
        return self._ref_structure

    @property
    def _is_rotated(self):
        return not np.isclose(self.prism.R, np.eye(3)).all()

    @property
    def positions(self):
        return self.from_voigt(self.ref_structure.positions, degree=1)

    @positions.setter
    def positions(self, new_positions):
        self.ref_structure.positions = self.to_voigt(new_positions, degree=1)

    @property
    def cell(self):
        return self.prism.A.flatten()[[0, 4, 8, 7, 6, 3]]

    @cell.setter
    def cell(self, new_cell):
        self.ref_structure.cell = self._prism.unfold_cell(new_cell)

    def from_voigt(self, values, degree=1):
        if degree == 1:
            if not self._is_rotated:
                return values
            return np.einsum('ij,...j->...i', self.prism.R, values)
        elif degree == 2:
            v = np.asarray(values)
            arr = [
                v[..., 0], v[..., 3], v[..., 4],
                v[..., 3], v[..., 1], v[..., 5],
                v[..., 4], v[..., 5], v[..., 2]
            ]
            if not self._is_rotated:
                return arr
            return np.einsum(
                'iI,jJ,...IJ->...ij',
                self.prism.R,
                self.prism.R,
                np.stack(arr, axis=-1).reshape(v.shape[:-1] + (3, 3))
            )
        else:
            raise ValueError('Only degree 1 or 2 is supported')

    def to_voigt(self, values, degree=1):
        if degree == 1:
            if not self._is_rotated:
                return values
            return np.einsum('ji,...j->...i', self.prism.R, values)
        elif degree == 2:
            # Symmetrization - mostly redundant
            v = 0.5 * (np.einsum('...ij->...ji', values) + values)
            if self._is_rotated:
                v = np.einsum(
                    'Ii,Jj,...IJ->...ij',
                    self.prism.R,
                    self.prism.R,
                    v
                )
            return v.reshape(v.shape[:-2], 9).T[[0, 4, 8, 7, 6, 3]].T
        else:
            raise ValueError('Only degree 1 or 2 is supported')

(to be honest I should have created a smaller example first...)

With this, in LAMMPS we completely stop using self.structure and instead all communication goes via self.structure._voigt, e.g. getting/setting positions (via self.structure._voigt.positions = new_positions). I understand that this doesn't fully solve the problem but at least this allows us to get rid of a lot of confusions.

[samwaseda]

Here a bit of explanation: structure._voigt is a shadow class of structure, which gives cell in the pseudo-Voigt notation (i.e. lxx lyy lzz lyz lxz lxy) and the positions accordingly. The setters allow for a translation from positions in the pseudo-Voigt notation to the normal ones.