Inconsistent energy and forces on conversion from ASEAtoms to Atoms

[jameskermode]

Consider the script

using JuLIP
import ASE

at0 = ASE.read("in1.xyz")
at1 = Atoms(at0)

set_constraint!(at0, FixedCell(at0))
set_constraint!(at1, FixedCell(at1))

sw = StillingerWeber()

@show norm(positions(at0) - positions(at1))
@show norm(cell(at0) - cell(at1))

@show energy(sw, at0) - energy(sw, at1)
@show norm(forces(sw, at0) - forces(sw, at1))

Output is

norm(positions(at0) - positions(at1)) = 0.0
norm(cell(at0) - cell(at1)) = 0.0
energy(sw, at0) - energy(sw, at1) = 1.0428900716501062
norm(forces(sw, at0) - forces(sw, at1)) = 9.507177398733772

i.e. positions and cell match but not computed quantities. The forces and energy of the ASEAtoms are the correct ones, those on the converted Atoms are wrong.

Am I doing something incorrect or is this a bug?

[jameskermode]

Forgot the input file:

$ cat in1.xyz
8
Lattice="4.85698085028444 1.15271546006108 -0.255327365498783 -1.41070126834903 3.42365804514259 0.112948147821099 -0.825925121792885 -1.79997392135745 3.84891948811008" Properties=species:S:1:pos:R:3:initial_magmoms:R:1:Z:I:1:castep_labels:S:1 pbc="T T T"
C       0.19575252       2.89628800       0.69082030       0.00000000        6 NULL
C      -0.74890209       2.87950453       1.81084364       0.00000000        6 NULL
C       1.87850136       0.89300291       0.40074928       0.00000000        6 NULL
C       2.83768421       4.21357737       1.27474085       0.00000000        6 NULL
C       2.94111961       0.96590247       1.56547879       0.00000000        6 NULL
C       1.23756824      -0.00377146       3.02548870       0.00000000        6 NULL
C       2.43824232       0.74214909       2.79634463       0.00000000        6 NULL
C       1.01883700       1.98038578       0.20645812       0.00000000        6 NULL

[cortner]

Lets assume it is a bug. How urgent?

[jameskermode]

Not terribly - working fine if I stick with ASEAtoms, by didn't want it to get forgotten.

[jameskermode]

One more datapoint: I checked pbc(at0) == pbc(at1), so I think it's the neighbour list which differs. Reason I think the ASEAtoms is correct is that minimisation converges, while the Atoms gets stuck.

[jameskermode]

It's a non-orthorhombic unit cell, so bugs in neighbourlist are probably to be expected.

[cortner]

I debugged the new nlist against matscipy general cell shapes, so I’m a bit surprised. I’ll have a look

[cortner]

With

at1 = rattle!(bulk(:C) * 2, 0.02)
at0 = ASE.ASEAtoms(at1)

the results are identical, so I'm not sure it is the shape of the cell

[cortner]

same if I randomly perturb the cell shape

[cortner]

maybe I have somewhere implicitly assumed that all positions are inside the cell? (which I think fails for your file)

(just talking out loud . . .)

[cortner]

yep - that seems to be it, I'll investigate further. . .

[cortner]

using JuLIP
import ASE

at0 = ASE.read("in1.xyz")
at1 = Atoms(at0)

sw = StillingerWeber()

@show norm(positions(at0) - positions(at1))
@show norm(cell(at0) - cell(at1))
@show pbc(at0) == pbc(at1)

@show energy(sw, at0) - energy(sw, at1) # 1.04289

X = positions(at1)
C = cell(at1)
X[2] -= C[2,:]
X[4] -= C[2,:]
set_positions!(at1, X)
@show energy(sw, at0) - energy(sw, at1) # 5.684341886080802e-14

[cortner]

I am drawing a blank: the rows of cell are the cell vectors, correct? This means that to get the barycentric coordinates I need y = inv(C)' * x?

[jameskermode]

Not sure quite what you mean by barycentric - it should be t = inv(C') * x for the fractional coordinates - transpose first, then inverse second.

[cortner]

That should actually be the same. Thanks.

[cortner]

I pushed a workaround to master; basically I'm just wrapping all atoms into the cell. But this is technically not fixing the bug.

@jameskermode when you have a minute could you check though whether this works for you now?

[cortner]

I am closing this since it is a NeighbourLists.jl bug: NeighbourLists.jl:#6

[jameskermode]

Thanks - I don't see the new commit yet, should it be in JuLIP or NeighbourLists master?

[cortner]

Neighbourlists.jl

[cortner]

@jameskermode I’ve publish new versions of JuLIP and NeighbourLists