Not orthorhombic unit cell

[hlscalon]

I am porting some code from python (mostly using ASE) to this library, and it's been great.

But I am stuck in a problem I cant seem to resolve.
When working with a FCC lattice with ASE, using a 3x3 vector like:

[ [5.1053 0.0 0.0] [0.0 8.8426 0.0] [0.0 0.0 12.5054] ]

I can port it to Aboria searching functions using min, max like:

vdouble3 min = vdouble3(0.0, 0.0, 0.0);
vdouble3 max = vdouble3(5.1053, 8.8426, 12.5054);
vbool3 periodic = vbool3(true, true, false);

particles.init_neighbour_search(min, max, periodic);

And it works great, result is the same as in ASE (using functions like get_distances, etc.)

But, when working with a HCP lattice (ie. 90º 90º 120º), like:

[ [8.8500 0.0 0.0] [-4.4249 7.6643 0.0] [0.0 0.0 9.3692] ]

Making stuff like:

vdouble3 min = vdouble3(-4.4249, 0.0, 0.0);
vdouble3 max = vdouble3(8.8500, 7.6643, 9.3692);

or

vdouble3 min = vdouble3(-4.4249, 0.0, 0.0);
vdouble3 max = vdouble3(8.8500, 8.8500, 9.3692);

I can't get the same result as in python.
Is it possible to use such unit cells (ie. not 90º 90º 90º) in Aboria ?

I am sorry if I am missing something trivial, not my area

[martinjrobins]

No, you didn't miss anything. Unfortunately Aboria is currently only able to use a hypercube domain (i.e. 90º 90º 90º), simply because I have never needed anything else. I would be open to enabling different types of domain, but I'm not sure how much I would need to adapt the existing spatial data structures to use, for example, a HCP lattice. My background is not in traditional molecular dynamics codes like ASE, so I'm not sure how this is normally implemented. I'll have a look around, but if you have any resources/links I can look at this would be helpful.

[hlscalon]

Thanks for your quick response!

Not sure if that helps, but the code in ASE that I think it's more closely related with this can be found in (taking for example the distance between the atoms):

https://wiki.fysik.dtu.dk/ase/_modules/ase/atoms.html#Atoms.get_distances
https://wiki.fysik.dtu.dk/ase/_modules/ase/geometry/geometry.html#get_distances
https://wiki.fysik.dtu.dk/ase/_modules/ase/geometry/geometry.html#find_mic
https://wiki.fysik.dtu.dk/ase/_modules/ase/geometry/cell.html#complete_cell

I can provide a MVCE on python (ASE) and Aboria (with fcc and hcp), if that can help you.

[martinjrobins]

Thanks @hlscalon , I'll have a look through these and get back to you (some time next week probably)

[martinjrobins]

Please do add the MVCE, that would be helpful

[hlscalon]

I created a simple gist:

https://gist.github.com/hlscalon/ccd0ce084532973d25c8835baffa19a7

I've hardcoded the min, max part in the cpp file. Then you can call with a random point:

./nnsearch fcc.xyz 1.5 1.5 1.5
python2 nnsearch.py fcc.xyz 1.5 1.5 1.5

[martinjrobins]

ok, so as far as I can tell these different lattices are just setting up a skew coordinate system with a transformation matrix A, so that

A * u = x

where u is the un-transformed vector position, and x is the transformed vector position (i.e. in the skew coordinate system)

my proposed solution to this is to assume that all the positions in Aboria's particle set are the un-transformed coordinates u, so if a user wants to get the real transformed coordinates x they simply multiply by the matrix A

The complication comes in with the neighbour search, where Aboria tries to search in the un-transformed space. So I would propose that Aboria's searches optionally take a Transform argument that is a function objects that the search can use internally to transform the coordinates to the correct, transformed space.

So the code to do a neighbour search would be something like

auto mult_A = [](vdouble3& x) { \\ set x = A * x }
auto mult_invA = [](vdouble3& x) { \\ set x = A^{-1} * x }

for (int i; i < particles.size(); ++i) {
     get<position>(particles)[i] = position_in_transformed_space
     mult_invA(get<position>(particles)[i]); // this changes the position to the un-transformed space
}

auto min = vdouble3(0.0, 0.0, 0.0);
auto max = vdouble3(1.0, 1.0, 1.0); \\ use the min/max coordinates of the un-transformed space
auto periodic = vbool3(true, true, false);

particles.init_neighbour_search(min, max, periodic);

// pass euclidean_search each the transformation lambda so it can use it internally
for (auto i = euclidean_search(particles.get_query(), rp, radius, mult_A); i != false; ++i) }
    // i.dx() here will return the vector between the two particles in the *transformed* space
    std::cout << "Found a particle with dx = " << i.dx() << " and id = " << get<id>(*i) << "\n";
}

what do you think @hlscalon, would this be suitable for what you need?

[hlscalon]

Thanks @martinjrobins for looking into this. I'll make some more tests but it seems fine for me.

[martinjrobins]

Great, I'll finish implementing this then and mark this issue when its ready for testing.