JuLIP.jl
JuLIP.jl copied to clipboard
JuliaMolSim
`read_extxyz` does something strange with the virials for 3-atom cells
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]])
@jameskermode -- any idea where this comes from? Should be quick to fix, but I wanted to check with you first.
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.
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.
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?
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.
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 https://github.com/JuliaMolSim/JuLIP.jl/issues/166#issuecomment-1636702132, or unsubscribe https://github.com/notifications/unsubscribe-auth/ACXHHTVRAMJVSUTPASFGYHDXQJD27ANCNFSM6AAAAAA2KVNJVQ . You are receiving this because you authored the thread.Message ID: @.***>