NaN energy when calling Atomistica potentials via ASE

[jameskermode]

The following script demonstrates the error. We tracked it down as far as set_position_dofs!().

atoms_broken_relaxed.xyz.txt

using PyCall
using JuLIP
@pyimport atomistica
@pyimport ase

atoms = ASEAtoms(ASE.ase_io.read("atoms_broken_relaxed.xyz"))
calc = ASECalculator(atomistica.TersoffScr())
set_calculator!(atoms, calc)

set_constraint!(atoms, FixedCell(atoms))
println(energy(atoms))
println(energy(atoms, dofs(atoms)))

Results are:

julia> include("run.jl")
-2040.7844479506173
NaN

i.e. second call to energy(atoms, dofs(atoms)) gives a NaN result.

[jameskermode]

Here's a slightly simpler but equivalent example:

using PyCall
using JuLIP
@pyimport atomistica
@pyimport ase

atoms = ASEAtoms(ASE.ase_io.read("atoms_broken_relaxed.xyz"))
calc = ASECalculator(atomistica.TersoffScr())
set_calculator!(atoms, calc)
set_constraint!(atoms, FixedCell(atoms))

println(energy(atoms))
set_positions!(atoms, positions(atoms))
println(energy(atoms))

[jameskermode]

Changing the definition of positions(at::ASEAtoms) from

positions(at::ASEAtoms) = copy( pyarrayref(at.po["positions"]) ) |> vecs

to

positions(at::ASEAtoms) = at.po[:get_positions]()' |> vecs

appears to solve this problem. I don't fully understand why, but since a copy is made in either case, avoiding the unsafe_wrap() is probably good. If there are no objections I will check this in after a bit more testing.

[cortner]

The fix makes 2 copies. Ok to make sure t works but we should find the bug.

[jameskermode]

Yes, I just realised that. I'll try with PyReverseDims().

[jameskermode]

I'm homing in on the problem. If a = PyArray(atoms.po["positions"]), then unsafe_wrap(Array, a.data, reverse(a.dims)) (what we use now) does not give array elements in the same order as transpose(unsafe_wrap(Array, a.data, a.dims)) (what we want) since the Python array is in C (row-major) order. I don't have a solution that doesn't involve a copy yet.

[jameskermode]

The plot thickens... I finally understand why this happens when reading from an .xyz file but not when creating structures with bulk(). In the former case the numpy positions array is Fortran-contiguous, in the latter it's C-contiguous. Perhaps a bug in the ASE xyz reader (it's me that wrote that).

[cortner]

very interesting.

[cortner]

There is this issue, which discusses a TransposeMatrix type, but I don't think it is implemented yet. But once it exists we can simply check what the storage type of the Python array is and import the unsafe_positions either as Matrix or as TransposeMatrix. Until then, maybe we need to drop unsafe_positions or make the behaviour undefined (as make it undefined whether one gets a reference or a copy).

[jameskermode]

I think it's just a bug in the ASE XYZ reader: I assemble the columns using vstack() which leads to Fortran memory ordering. I've created a pull request with a fix, once the ASE guys merge this then we can revert commit 66bcf2e and close this issue.

[cortner]

thats fine for me if there is a guarantee that ASE will always use the memory layout we assume now? It just never occurred to before that it might not, but now I am worried.

[jameskermode]

I think it should always be in C order, but the problem is that from Python you often don't notice if the memory is in C or Fortran order since numpy abstracts it completely, unless you look very carefully (or notice slow downs).

[jameskermode]

If we are worried we could actually check the flags of the numpy array before doing the unsafe_wrap. I'll have a go at this.

[cortner]

this is my point - can we guarantee that ASE will always have positions stored in C order? If not then we need to check before we import them into JuLIP.

[cortner]

right, that's what I meant. And maybe for the time being, if positions are stored in F order then let unsafe_ thing just return a copy instead of a reference?

[jameskermode]

Probably unsafe_positions should raise an error rather than returning a copy, otherwise you might be expecting to be able to modify in place. But positions() should be made to return a copy in correct memory layout.

[cortner]

Ok then let's raise an error for now and switch it to TransposeMatrix once it exists. C

…

Sent from my iPhone On 16 Feb 2017, at 21:16, James Kermode <notifications@github.com<mailto:notifications@github.com>> wrote: Probably unsafe_positions should raise an error rather than returning a copy, otherwise you might be expecting to be able to modify in place. But positions() should be made to return a copy in correct memory layout. — You are receiving this because you commented. Reply to this email directly, view it on GitHub<#48 (comment)>, or mute the thread<https://github.com/notifications/unsubscribe-auth/AHVRQwohJLes_AGbzQnlK1gVSOCRiLF-ks5rdLy6gaJpZM4MDSZe>.

[cortner]

I guess this will be fixed once we create a pure Julia Atoms.

[cortner]

closing for now - we can reopen if it comes up again