`read_extxyz` does something strange with the virials for 3-atom cells

[wcwitt]

For 3-atom cells, read_extxyz converts the virials to a vector of vectors. For all other cell sizes, read_extxyz leaves the virials as a matrix. This behaviour was causing some weird things in ACE1pack.

I believe this function explains it: https://github.com/JuliaMolSim/JuLIP.jl/blob/3d23cd8a4a30d46643300da8821c0ee9f953010a/src/atoms.jl#L419

Minimal example: consider this train.xyz file.

5
Lattice="1.92847 0.0 0.0 -0.134535 5.944528 0.0 -0.869345 -2.444608 6.542313" Properties=species:S:1:pos:R:3:forces:R:3 energy=-798.7837161363719 virial="64.28957258538264 1.273211496227935 0.7897628948211044 1.273211496227935 41.18384738457577 -5.666667903558281 0.7897628948211044 -5.666667903558281 4.8516786004571895" pbc="T T T"
Si       0.60663610      -1.09674126       6.08761291      -0.24701513      -0.49480659      -0.58493040
Si      -0.01657097      -0.27774239       3.73073884      -0.23107425     -10.72330762       0.81933772
Si       0.04032046       1.41472773       3.62884558      -0.18807675       6.12381475       0.63003909
Si       0.32032607       3.37283612       3.42038853       1.03007300       3.42957453       2.90447579
Si      -0.20521643       4.49652268       1.71143436      -0.36390687       1.66472493      -3.76892221
3
Lattice="3.38243 0.0 0.0 -1.34636 3.536453 0.0 -0.567482 -0.944697 3.761965" Properties=species:S:1:pos:R:3:forces:R:3 energy=-482.0028593517148 virial="10.944787239437682 9.4970514154512 -8.623491235978445 9.4970514154512 12.910487855807412 -0.381666084697965 -8.623491235978445 -0.381666084697965 19.641593532042016" pbc="T T T"
Si      -0.17616872      -0.35226186       2.87414811      -3.25470768      -2.68899422      -2.43446614
Si       0.71333136       1.97408966       0.76444673      -6.55425455      -0.30111898      12.95430646
Si       1.02119970       1.19876219       3.07014725       9.80896223       2.99011320     -10.51984031
4
Lattice="2.88569 0.0 0.0 -0.638994 3.552377 0.0 -0.112749 -1.442029 5.853055" Properties=species:S:1:pos:R:3:forces:R:3 energy=-621.183963984326 virial="26.12667767531235 -14.556200330474802 3.3416789123311847 -14.556200330474802 20.374956390145357 10.082086344818723 3.3416789123311847 10.082086344818723 8.225262348140962" pbc="T T T"
Si       0.60890556       0.72128915       1.99135669      13.35110884     -19.76196861      -7.41608320
Si       2.34457095       1.59319193       1.90416841     -15.59357600       6.27928433      -7.68632778
Si       2.14974448      -1.05746357       5.68069857      -0.26542886      -0.46148720       0.44900284
Si       0.08814170       1.87266310       2.74974751       2.50789602      13.94417148      14.65340814
3
Lattice="2.99718 0.0 0.0 -0.343807 3.212947 0.0 -1.293007 -1.371618 4.673009" Properties=species:S:1:pos:R:3:forces:R:3 energy=-484.41007978749445 virial="8.770665763052879 8.460719982108959 2.711326048251247 8.460719982108959 15.302936481334566 3.589458205403738 2.711326048251247 3.589458205403738 6.950102933550505" pbc="T T T"
Si       1.80763774       1.59699271       2.82646784      -0.27266431      -0.12002321       1.79161286
Si       0.26785904       2.74583519       1.14218190       6.66691665       9.03722027       3.14394393
Si       2.25383832       1.33428485       0.64319617      -6.39425233      -8.91719706      -4.93555679

The following

using JuLIP

train = read_extxyz("train.xyz")
for t in train
    @show t.data["virial"]
end

produces

t.data["virial"] = JuLIP.JData{Float64}(Inf, 0.0, [64.28957258538264 1.273211496227935 0.7897628948211044; 1.273211496227935 41.18384738457577 -5.666667903558281; 0.7897628948211044 -5.666667903558281 4.8516786004571895])
t.data["virial"] = JuLIP.JData{Float64}(Inf, 0.0, StaticArraysCore.SVector{3, Float64}[[10.944787239437682, 9.4970514154512, -8.623491235978445], [9.4970514154512, 12.910487855807412, -0.381666084697965], [-8.623491235978445, -0.381666084697965, 19.641593532042016]])
t.data["virial"] = JuLIP.JData{Float64}(Inf, 0.0, [26.12667767531235 -14.556200330474802 3.3416789123311847; -14.556200330474802 20.374956390145357 10.082086344818723; 3.3416789123311847 10.082086344818723 8.225262348140962])
t.data["virial"] = JuLIP.JData{Float64}(Inf, 0.0, StaticArraysCore.SVector{3, Float64}[[8.770665763052879, 8.460719982108959, 2.711326048251247], [8.460719982108959, 15.302936481334566, 3.589458205403738], [2.711326048251247, 3.589458205403738, 6.950102933550505]])