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.

[martinjrobins]

Hi @hlscalon. This work is now finished and merged with the master branch. Are you able to have a look at the documentation link below and test this out on your problem?

https://martinjrobins.github.io/Aboria/aboria/neighbourhood_searching.html#aboria.neighbourhood_searching.coordinate_transforms

[hlscalon]

Hi @martinjrobins, thanks for the time taking looking into this. I have some doubts though.

Let's say I have the file bcc.xyz (attached MVCE.zip).

It's unit cell is

[ 5.76  0.00  0.00 ]
[ 0.00  5.76  0.00 ]
[ 4.32  4.32  4.32 ]

So it's transformation matrices would be this, right ?

auto mult_A = [](const vdouble3 & v) {
    return vdouble3(v[0] + 4.32 * v[2], v[1] + 4.32 * v[2], v[2]);
};

auto mult_invA = [](const vdouble3 & v) {
    return vdouble3(v[0] - 4.32 * v[2], v[1] - 4.32 * v[2], v[2]);
};

But what about the min/max vectors ? Early in this gist I have used for the 'fcc' file the following values (5.10.. is not the max x point in that file, for example, which is in fact 3.82.. And that is ok, I am using several other files with that structure correctly):

vdouble3 min = vdouble3(0.0, 0.0, 0.0);
vdouble3 max = vdouble3(5.105310960166873, 8.842657971447274, 12.505406830647296);

But now, if I try to use the 'un-transformed' value of [5.76, 5.76, 4.32] for max, for example, i got [-12.9024, -12.9024, 4.32], which is unlikely to be correct. For testing I am using the same gist as before, just modified for using the skew coordinations (transforming the points of the file with mult_invA and then using mult_A for to transform it back. File attached MVCE.zip)

(I am probably missing something very trivial, but I couldn't realize what it is)

[martinjrobins]

What is the definition of these numbers?:

[ 5.76 0.00 0.00 ] [ 0.00 5.76 0.00 ] [ 4.32 4.32 4.32 ]

To use the transforms in Aboria, you will need a matrix that maps vectors in the untransformed space (i.e. a unit cube from (0,0,0) to (1,1,1)) to your new skewed space. So if A is your transform matrix, then A*(1,1,1)^T should give the maximum extent of your skewed space. But it looks like just treating the numbers you give above will not do this?

[martinjrobins]

Found this set of slides which might be useful:

http://cloud.crm2.univ-lorraine.fr/pdf/Antwerp2016/Aroyo.pdf

So if you have a look at the coordination transformation section in the slides, the matrix A that I referred to earlier is the matrix P in the slides