Coordinate transformations in lammps simulations with changing cells

[Leimeroth]

In lammps dump parsing coordinates and cell are transformed between lammps and pyiron coordinate systems using: ''' prism = self._prism rotation_lammps2orig = self._prism.R.T ''' As far as I can tell self._prism creates all rotation matrices etc. based on the initial cell. These matrices are than used to transform coordinates and cell shape between the coordinate systems. I assumed that this transformation has to be done by creating a prism and corresponding matrices based on the cell of the current step instead of the initial cell, especially when the cell shape changes.

[Leimeroth]

maybe @liamhuber or @samwaseda could explain what is going on there?

[jan-janssen]

I feel this issue is redundant https://github.com/pyiron/pyiron_atomistics/issues/673

[Leimeroth]

I feel this issue is redundant #673

673 is clearly a bug. In this case I simply do not understand whether this works as it is supposed to or if all lammps jobs with changing cell do something wrong during parsing by not updating the prism used for coordinate transformations for each step.

[liamhuber]

I assumed that this transformation has to be done by creating a prism and corresponding matrices based on the cell of the current step instead of the initial cell, especially when the cell shape changes.

I have very low confidence in my answer here and would recommend explicitly creating a scenario where such a transformation takes place and seeing if it still works. With that said, I think the initial prism should be fine. The prism just transforms coordinates from whatever pyiron is using to an upper(lower?) triangular bravais matrix. Once that's done initially, I would expect evolutions of the cell within the triangular lammps framework to continue reverse-transforming back to pyiron ok? I.e. if $cell{pyiron}[0] = R^{-1} * cell{lammps}[0]$, and $cell{lammps}[n] = cell{lammps}[0] + dcell{lammps}$, shouldn't we still just get $cell{pyiron}[n] = R^{-1} cell{lammps}[n] = R^{-1} (cell{lammps}[0] + dcell_{lammps})$

[samwaseda]

I agree with @liamhuber.

This being said, I want to remind you that we might have a situation, where the LAMMPS might give rotation-free shape changes, while through this transformation pyiron might create a rotation (because the x-axis is rigid, and the y-axis also for the z-axis). I don't think there's anything we could do about it, but we shouldn't forget this because I'm sure that one day someone is gonna have a problem and we don't understand where it comes from...

[Leimeroth]

cellpyiron[0]=R−1∗celllammps[0], and celllammps[n]=celllammps[0]+dcelllammps, shouldn't we still just get cellpyiron[n]=R−1celllammps[n]=R−1(celllammps[0]+dcelllammps)

Shouldn't R be R(n)? And is this what is actually done? I interpreted the code in dump parsing to calculate pyiron positions based on the relative positions as given in the dump file. These positions are relative to the current cell and I don't see where dcell comes in to the calculation to correct for the changes between cell[0] and cell[n]

[liamhuber]

Shouldn't R be R(n)?

Perhaps I misunderstood -- for interactive jobs, where the structure is being re-set, then I agree that the rotation matrix is step-dependent. Then I think this comes back to issue #673 that Jan mentioned.

In case it is just a regular lammps job in which the cell shape changes during the Lammps run, I still think we need only one rotation matrix:

• We start with a pyiron Bravais matrix $M_0$
• We construct a rotation matrix that makes a coordinate transform so that this matrix is triangular $M'_0 = RM_0$
• Lammps runs and the cell evolves, but of course stays triangular, $M'_n$
• We use the inverse of the original transformation to get from the triangular coordinate frame back to our original coordinate frame $M_n = R^{-1} M'_n$

$dcell$ was just my nomenclature for the matrix that tracks the difference between our original and final triangular matrices, $dcell := cell[n] - cell[0]$. My argument is that if $R^{-1}$ maps $M'$ to $M$, then it also maps $dM'$ to $dM$ just fine, since it $R$ was created to map between a particular full-matrix coordinate frame and a particular triangular coordinate frame. $M'_n$ and $M'_0$ are in the same coordinates, so I figure the coordinate map should stay valid. However, I have always found geometry surprisingly difficult and I am open to the possibility that this is not correct.

If we're talking about relative ($r$) and absolute ($p$) positions, then $p = Mr$, so $p_n = M_n r_n = R^{-1} M'_n r'_n = R^{-1} p'_n$ -- i.e. we will definitely need the final cell to map relative positions back to final positions, but the rotation matrix (prism) is static. I can't vouch that we use rotation_lammps2orig everywhere we ought to, but I'm comfortable with it not getting updated.