make unitcell matrix column major

[lmiq]

The docs (here) imply that the matrix returned by matrix is column major. However, this is not what is actually happening. (example file attached).

Example:

julia> traj = Chemfiles.Trajectory("./traj_comp.xtc");

julia> frame = Chemfiles.read_step(traj,0);

julia> unitcell = Chemfiles.UnitCell(frame);

julia> Chemfiles.matrix(unitcell)
3×3 Matrix{Float64}:
 77.395     0.0     0.0
 38.6975   67.026   0.0
 38.6975  -22.342  63.1928

as such, this matrix does not satisfy the conditions gromacs requires.

Just transposing the matrix would be easy to solve as Matrix(transpose(Chemfiles.matrix(...))). However, I'm not sure if it makes the definition inconsistent with the return of the unit cell in terms of lengths and angles. For instance, in this trajectory, we get:

julia> Chemfiles.matrix(unitcell)
3×3 Matrix{Float64}:
 77.395     0.0     0.0
 38.6975   67.026   0.0
 38.6975  -22.342  63.1928

julia> Chemfiles.lengths(unitcell)
3-element Vector{Float64}:
 77.39500045776367
 77.39500188456586
 77.39500082545612

julia> Chemfiles.angles(unitcell)
3-element Vector{Float64}:
 90.00000019466002
 60.000000157157025
 60.000000609835716

But if we tranpose the matrix, we do not recover these angles and lengths:

julia> mt = Matrix(transpose(Chemfiles.matrix(unitcell)));

julia> ut = Chemfiles.UnitCell(mt)
Chemfiles.UnitCell(Chemfiles.CxxPointer{Chemfiles.lib.CHFL_CELL}(Ptr{Chemfiles.lib.CHFL_CELL} @0x000000000910c100, false))

julia> Chemfiles.angles(ut)
3-element Vector{Float64}:
 108.4349482115009
  65.90515744788931
  75.03678223987454

julia> Chemfiles.lengths(ut)
3-element Vector{Float64}:
 94.78912988199575
 70.65164746673275
 63.19275379180908

The simplest alternative would be to fix the docs and indicate that the matrix is row-major. That is somewhat inconsistent with Julia, but then at least it lets to the user the decision to transpose the matrix and arch with the consequences.

The other alternative is to check if this PR is fine, I have only transposed the matrix back and forth when converting it from the UnitCell.

traj_comp.xtc.zip

[Luthaf]

Right, this is an issue that comes from the C++ side of things. Over there, the data is stored as row major, but containing the cell vectors as columns (which is confusing, and which I want to change in the next breaking release).

When passing data to Julia, the data is interpreted as column major, meaning the cell vectors are now rows of the matrix ...

---

I think this PR is a nice way of handling this disparity until we update the representation on the C++ side. I'm not sure if I should put these changes into a new release though, since users of Chemfiles.jl might already be handling the rows/column disparity, and this would be a breaking change for them. So you would need to pull in chemfiles from git to have these changes.

The alternative would be updating the docs for now, and waiting for the next breaking release to change the behavior. Do you have a preference about this?

[lmiq]

I thought about that. It is not clear to me if that can be considered a bug and the docs are correct, though.

I don't have a strong opinion. Feel free to deal with it in the way you find safer.

[lmiq]

I sort of need a definition on this. Lets keep it as is, and let the change in row order for the next breaking version? Then I'll handle that within my package.

[Luthaf]

That works for me. Sorry, I should have told you I'm on holidays right now and wanted to address this after I'm back.