`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:

JuLIP.jl/src/atoms.jl

Line 419 in 3d23cd8

_read_convert(value::AbstractMatrix, natoms::Int) = size(value, 2) == natoms ? vecs(value) : value

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]])

[cortner]

@jameskermode -- any idea where this comes from? Should be quick to fix, but I wanted to check with you first.

[jameskermode]

Yes, I understand the problem and Chuck’s link is correct. The problem is that JuLIP’s additional Atoms data (unlike ASE) does not distinguish between per-atom and per-config data, so I have to guess. The heuristic is that if an array is 3 * N_atoms it’s per atom vector data, otherwise it’s per-config. This fails for 3-atom configs. Not sure of the best fix, maybe treating virial as a special case since it will always be per-config is enough for this particular issue.

[jameskermode]

Actually, we do know the origin when reading, so could special case that _read_convert() function to work differently for info (per config) and arrays (per atom). Writing is more difficult once it’s all been thrown into one Atoms.data bucket.

[cortner]

that's an interesting problem. nasty. maybe in practice writing is ok since we treat the cell as a matrix.

But if we wanted to do this systematically we could add a label to atoms-data to specify whether it is per atom or per-config. That should be easy to do. Would you like me to start a PR for this?

[jameskermode]

Labelling the data would be the best solution. If you have time to make a PR that would be great, I’m happy to check it over.

[wcwitt]

Thanks, both. There’s now a virial-specific workaround in ACE1pack, so no urgency from my end. But I’ll happily test and replace it when you’re ready

…

On Sat, 15 Jul 2023 at 08:52, James Kermode ***@***.***> wrote: Labelling the data would be the best solution. If you have time to make a PR that would be great, I’m happy to check it over. — Reply to this email directly, view it on GitHub <#166 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ACXHHTVRAMJVSUTPASFGYHDXQJD27ANCNFSM6AAAAAA2KVNJVQ> . You are receiving this because you authored the thread.Message ID: ***@***.***>