Polycrystal orientations are not completely random

[jfikar]

Dear Pierre,

I think I found another problem in polycrystal. This time I was creating equilateral randomly oriented grains of tungsten and I was looking at its elastic properties.

box 50 50 50
random 12

I did 1000 of such samples and the Young modulus in Z is on average for some reason smaller than in X and Y. I could also get the standard deviation and the difference of Young modulus in Z and X or Y is 20x the standard deviation. So probably it is not a random effect.

Moreover, tungsten has quite isotropic elastic constants. I think the effect would be grater for anisotropic material.

As a test I created one big sample, equal to my 1000 smaller samples

box 500 500 500
random 12000

And I inspected the angles written in poly_param.txt, the first and third angle are uniformly distributed -180..180deg and the second 0..180deg with maximum in 90deg.

I think I found the place related to this problem in mode_polycrystal.f90, it is around line 775 you have a comment:

!NOTE: just multiplying the three random numbers by 2pi results in a biased distribution
!See J.J. Kuffner, Proc. 2004 IEEE Int'l Conf. on Robotics and Automation (ICRA 2004)
       !or http://mathworld.wolfram.com/SpherePointPicking.html

This is true, but it seems, you are mixing the first method (eq. 1 and 2 referenced by Wolfram), where just two angles are generated to cover the sphere (latitude ϕ and longitude θ) with a third independent rotation. I'm not sure it is correct.

I did check the rotation matrices (similar to the ones you have in the code) for the angles from poly_param.txt and plotted the position of X-axis (1 0 0) after the 12000 random rotations and I get only half of the sphere covered with a lot of points around the pole and not so much at the equator.

I think this is the reason of the observed anisotropy of the generated polycrystals. I think the proper random rotation is in fact more complex. Probably the best is the approach of Arvo. You'll get directly the rotation matrix. I'm however not sure, how to get the two angles back. One angle is obvious.

[pierrehirel]

Dear jfikar,

Thank you for pointing this out, and for your detailed report of the problem.

Indeed it appear that the implementation is wrong, however I think the idea outlined by Wolfram is correct. So, maybe the easiest way to fix this is to just set the third angle equal to zero? That would produce only two rotations but should cover the whole sphere. Can you try this, i.e. modify the following line:

randarray(o) = 0.d0

By the way, when running the code with debug option (add "-v 4" to your command-line), a file named "atomsk_angles.xyz" is produced. It contains the positions of the random points projected on a sphere, and after setting the third angle to zero, I find that the points are uniformly distributed on a sphere. Example with 12000 grains:

Regards

[jfikar]

Dear Pierre,

Wolfram is correct, however it solves a different problem. The formula in Wolfram is for randomly choosing just one vector, which will cover a sphere uniformly.

But we need a 3D rotation, we need a new basis of three vectors. Well, third one is of course easy to find once you have the two. Unfortunately, this problem is not possible to solve by first finding the direction covering a sphere uniformly and then add a random rotation.

I'm not sure, what exactly the points in atomsk_angles.xyz are. What I did is to construct the three rotation matrices from your angles in poly_param.txt, multiply them together, first applying the rotation around axis X, then Y and then Z: T=RzRyRx.

Then I take the old vertical basis vectors [1 0 0], [0 1 0], [0 0 1], rotate them all by T i.e. T*[1 0 0] etc. and look where the new basis will be projected on the sphere.

As it happens, Y and Z look close to correct, but X does not (that's the picture in my previous post).

Setting the third angle to zero makes the things worse. X looks now only like a half-circle. The first rotation along X does not change X. Second rotation is limited to 0-180 and third rotation is now set to zero, so it should be only a half-circle.

At the same time Y and Z endpoints after this rotation look good, sphere is almost uniformly covered. There is however the half-circle or circle missing and there are two denser points at ends of this half-circle.

Maybe if you look at your picture, or use a slightly less data points, you'll see the same.

Unfortunately I don't know about any easy method to solve this problem. Three rotations with random angles -180..180 are not correct either.

I have implemented the Arvo matrices in Matlab and this works well.

There is another possibility to do it by QR decomposition of normally (not uniformly) distributed 3x3 matrix. This works as well, but I needed three extra random numbers to flip randomly sign of the three columns of T. Otherwise the resulting endpoints of the rotated base are on a half-spheres. I think it may depend on the method actually used for QR decomposition.

I'm not sure if you need only the full rotation matrix (this I have), or you would like to have the three independent rotation angles. As I understood, that might be a problem.

Any 3D rotation matrix should have eigenvalues 1, and two complex conjugate. The eigenvector corresponding to the 1 is the rotation axis and the two complex eigenvalues give the rotation angle around this axis. Probably this is our third angle and first two are used to get to the eigenvector.

[pierrehirel]

OK, I think I understand the issue better now. Unfortunately I do not have much time at the moment to implement a correct algorithm in the code.

I guess the matrices are sufficient for the purpose of constructing the polycrystal. The angles are not necessary for the construction, but Atomsk writes them into a text file ("*_param.txt") so that the user knows them, and can reproduce the same polycrystal at a later time. So, it would be nice to have them.

Another solution could be to read rotation matrices instead of angles. I guess this could also work.

Anyway, if you have an implementation of Arvo method (or any working method), you may send it to me (by email or otherwise), or if you are willing, modify the code and ask for a pull request.

Thank you for your time dealing with this issue.

[jfikar]

I'm not completely sure, where exactly in your mode_polycrystal.f90 does the rotation take place. It seems to be in multiple locations, right?

However, the working Arvo code in Matlab looks quite simple:

random_num=rand(1,3)

theta=random_num(1)*2*pi;
phi=random_num(2)*2*pi;
z=random_num(3);
r=sqrt(z);

V=[r*cos(theta), r*sin(theta), sqrt(1-z)];
R=[cos(phi), sin(phi), 0; -sin(phi), cos(phi), 0; 0, 0, 1];
T=(2*V'*V-eye(3))*R

where random_num are three uniformly distributed random numbers 0..1 and T is the final 3D rotational matrix to be multiplied by column vectors like this: rotated_vec=T*old_vec

I can also provide 4 test cases:

random_num =

     0     0     0

T =

    -1     0     0
     0    -1     0
     0     0     1

random_num =

                       0.5                       0.5                       0.5

T =

     -2.22044604925031e-16      1.22464679914735e-16                        -1
      2.44929359829471e-16                         1      1.22464679914735e-16
                         1     -2.44929359829471e-16      2.22044604925031e-16

random_num =

     1     1     1

T =

                         1     -7.34788079488412e-16                         0
     -7.34788079488412e-16                        -1                         0
                         0                         0                        -1

random_num =

                       0.1                       0.2                       0.3

T =

        -0.459016994374947        -0.489403985718866         0.741476323045772
         0.842075137094349         0.026393202250021         0.538714082201795
        -0.283218753550188         0.871657695220709                       0.4

[pierrehirel]

Dear jfikar,

I finally found the time to implement the method from Arvo into the mode "--polycrystal" of Atomsk. It is now part of the master branch available from Github.

I tested it briefly and it seems to produce the desired distribution of random angles. When you have time, I would appreciate if you could also verify the results, and have a look at the code to validate the implementation.

Best regards

[jfikar]

Dear Pierre,

I've tried the new version and it looks good. The original code with 1000 samples with 12 grains of tungsten 50x50x50 gives average Young modulus: Yx=189.07±0.43, Yy=189.14±0.43, Yz=181.53±0.43 GPa

Note lower value of Yz. While the new code for the same test gives Yx=193.41±0.42, Yy=193.48±0.43, Yz=193.49±0.42 GPa

I've also tried the same with FCC copper, which is more anisotropic. The result for old code and 300 samples: Yx=69.94±0.42, Yy=69.78±0.43, Yz=65.86±0.43 GPa

Again Yz has lower value. And the new code: Yx=72.32±0.42, Yy=72.80±0.45, Yz=71.92±0.45 GPa.

So the new code works, thank you very much. I'll also look at the code later.

[pierrehirel]

Thank you for performing these tests. I leave this issue open for a time, so that you may comment on the code as well.

[jfikar]

As far as I can tell, your code on lines 804-825 is identical to my code in Matlab. I'm surprised the conversion from the rotation matrix back to three angles in so easy (your lines 916-918). But it also seems to work well.

I've just found one unstability point going from three random numbers to rotation matrix, then to three angles and construct the rotation matrix again from the angles (random_num → T → P → T2) for:

random_num =

                       0.5                       0.5                       0.5
T =

     -2.22044604925031e-16      1.22464679914735e-16                        -1
      2.44929359829471e-16                         1      1.22464679914735e-16
                         1     -2.44929359829471e-16      2.22044604925031e-16

P =
                   -47.8056204523404                       -90           132.19437954766

T2 =

     -1.90262369842273e-16        -0.995208235961456       -0.0977781523372533
      2.09871527663731e-16        0.0977781523372533        -0.995208235961456
                         1      -2.0987152766373e-16      1.90262369842273e-16

but changing any of the random numbers very slightly from 0.5 (+ or - 1e-9) gives again the same matrices T and T2. So it should not be a problem.

I also tried to visualise the endpoints of basal vectors after the rotation from the three angles written in poly_param.txt, the same as I did in the first post, and all of them seem fine i.e. covering sphere uniformly.

Good job and thank you