How to calculate quaternions from principal axes of an ellipsoid?

The name of the pictureThe name of the pictureThe name of the pictureClash Royale CLAN TAG#URR8PPP











up vote
1
down vote

favorite












I am working a project that involves calculating quaternions of ellipsoids providing that I already know the unit vectors of their principal axes. The reason for this is the program I am using to work with these ellipsoids require quaternions. The description from the manual of that program of quaternions is: "The particles must define a quaternion for their orientation. Note that particles defined as ellipsoid have 3 shape parameters. The 3 values must be non-zero for each particle set by this command. They are used to specify the aspect ratios of an ellipsoidal particle, which is oriented by default with its x-axis along the simulation box’s x-axis, and similarly for y and z. If this body is rotated (via the right-hand rule) by an angle theta around a unit rotation vector (a,b,c), then the quaternion that represents its new orientation is given by $(cos(theta/2), a*sin(theta/2), b*sin(theta/2), c*sin(theta/2))$"



From some source in the literature (which I am not sure if is correct), I calculate quaternions like this:



  • Construct a orientation matrix whose columns are the principal axes

  • Solve for e-values and e-vectors of that matrix. The e-values should be (1, $cos(theta) -i sin(theta)$,$cos(theta) +i sin(theta)$) where theta is the rotational angle. The unit e-vector $e = [e_1 e_2 e_3]^T$ that correspond to e-value 1 is an invariant principal axis of rotation.

  • The quaternions are:
    $q_w = cos(theta/2)$, $q_i = e_1sin(theta/2)$, $q_j = e_2sin(theta/2)$, $q_k = e_3sin(theta/2)$

However, when I try to check if the ellipsoids using the quaternions I calculated agree with the principal axes (using a graphic tool), they don't agree with each other.



I am not an expert in Maths, and only do this because of my project. I really appreciate if someone can help me to point out what I did wrong, or give me some good references, or even just explain to me what the manual description means.



Cheers,







share|cite|improve this question
























    up vote
    1
    down vote

    favorite












    I am working a project that involves calculating quaternions of ellipsoids providing that I already know the unit vectors of their principal axes. The reason for this is the program I am using to work with these ellipsoids require quaternions. The description from the manual of that program of quaternions is: "The particles must define a quaternion for their orientation. Note that particles defined as ellipsoid have 3 shape parameters. The 3 values must be non-zero for each particle set by this command. They are used to specify the aspect ratios of an ellipsoidal particle, which is oriented by default with its x-axis along the simulation box’s x-axis, and similarly for y and z. If this body is rotated (via the right-hand rule) by an angle theta around a unit rotation vector (a,b,c), then the quaternion that represents its new orientation is given by $(cos(theta/2), a*sin(theta/2), b*sin(theta/2), c*sin(theta/2))$"



    From some source in the literature (which I am not sure if is correct), I calculate quaternions like this:



    • Construct a orientation matrix whose columns are the principal axes

    • Solve for e-values and e-vectors of that matrix. The e-values should be (1, $cos(theta) -i sin(theta)$,$cos(theta) +i sin(theta)$) where theta is the rotational angle. The unit e-vector $e = [e_1 e_2 e_3]^T$ that correspond to e-value 1 is an invariant principal axis of rotation.

    • The quaternions are:
      $q_w = cos(theta/2)$, $q_i = e_1sin(theta/2)$, $q_j = e_2sin(theta/2)$, $q_k = e_3sin(theta/2)$

    However, when I try to check if the ellipsoids using the quaternions I calculated agree with the principal axes (using a graphic tool), they don't agree with each other.



    I am not an expert in Maths, and only do this because of my project. I really appreciate if someone can help me to point out what I did wrong, or give me some good references, or even just explain to me what the manual description means.



    Cheers,







    share|cite|improve this question






















      up vote
      1
      down vote

      favorite









      up vote
      1
      down vote

      favorite











      I am working a project that involves calculating quaternions of ellipsoids providing that I already know the unit vectors of their principal axes. The reason for this is the program I am using to work with these ellipsoids require quaternions. The description from the manual of that program of quaternions is: "The particles must define a quaternion for their orientation. Note that particles defined as ellipsoid have 3 shape parameters. The 3 values must be non-zero for each particle set by this command. They are used to specify the aspect ratios of an ellipsoidal particle, which is oriented by default with its x-axis along the simulation box’s x-axis, and similarly for y and z. If this body is rotated (via the right-hand rule) by an angle theta around a unit rotation vector (a,b,c), then the quaternion that represents its new orientation is given by $(cos(theta/2), a*sin(theta/2), b*sin(theta/2), c*sin(theta/2))$"



      From some source in the literature (which I am not sure if is correct), I calculate quaternions like this:



      • Construct a orientation matrix whose columns are the principal axes

      • Solve for e-values and e-vectors of that matrix. The e-values should be (1, $cos(theta) -i sin(theta)$,$cos(theta) +i sin(theta)$) where theta is the rotational angle. The unit e-vector $e = [e_1 e_2 e_3]^T$ that correspond to e-value 1 is an invariant principal axis of rotation.

      • The quaternions are:
        $q_w = cos(theta/2)$, $q_i = e_1sin(theta/2)$, $q_j = e_2sin(theta/2)$, $q_k = e_3sin(theta/2)$

      However, when I try to check if the ellipsoids using the quaternions I calculated agree with the principal axes (using a graphic tool), they don't agree with each other.



      I am not an expert in Maths, and only do this because of my project. I really appreciate if someone can help me to point out what I did wrong, or give me some good references, or even just explain to me what the manual description means.



      Cheers,







      share|cite|improve this question












      I am working a project that involves calculating quaternions of ellipsoids providing that I already know the unit vectors of their principal axes. The reason for this is the program I am using to work with these ellipsoids require quaternions. The description from the manual of that program of quaternions is: "The particles must define a quaternion for their orientation. Note that particles defined as ellipsoid have 3 shape parameters. The 3 values must be non-zero for each particle set by this command. They are used to specify the aspect ratios of an ellipsoidal particle, which is oriented by default with its x-axis along the simulation box’s x-axis, and similarly for y and z. If this body is rotated (via the right-hand rule) by an angle theta around a unit rotation vector (a,b,c), then the quaternion that represents its new orientation is given by $(cos(theta/2), a*sin(theta/2), b*sin(theta/2), c*sin(theta/2))$"



      From some source in the literature (which I am not sure if is correct), I calculate quaternions like this:



      • Construct a orientation matrix whose columns are the principal axes

      • Solve for e-values and e-vectors of that matrix. The e-values should be (1, $cos(theta) -i sin(theta)$,$cos(theta) +i sin(theta)$) where theta is the rotational angle. The unit e-vector $e = [e_1 e_2 e_3]^T$ that correspond to e-value 1 is an invariant principal axis of rotation.

      • The quaternions are:
        $q_w = cos(theta/2)$, $q_i = e_1sin(theta/2)$, $q_j = e_2sin(theta/2)$, $q_k = e_3sin(theta/2)$

      However, when I try to check if the ellipsoids using the quaternions I calculated agree with the principal axes (using a graphic tool), they don't agree with each other.



      I am not an expert in Maths, and only do this because of my project. I really appreciate if someone can help me to point out what I did wrong, or give me some good references, or even just explain to me what the manual description means.



      Cheers,









      share|cite|improve this question











      share|cite|improve this question




      share|cite|improve this question










      asked Aug 21 at 10:23









      Nguyen Thi Lan Huong

      82




      82




















          2 Answers
          2






          active

          oldest

          votes

















          up vote
          1
          down vote



          accepted










          Knowing the orientation vectors is essentially the same as knowing the transformation matrix, which is constructed by putting these vectors (they have to be normalized and orthogonal) into matrix columns.



          First possible issue: arrange the vectors so that they form a positively oriented triad (the determinant of the matrix should be 1, not -1). Otherwise, your matrix is not a rotation but a rotation+reflection. Fixing it is easy, as it's an ellipsoid and you can flip the sign of any main axis vector without changing anything. Second possible issue is that the order of axes is different and you are putting the long axis where the short axis should be (it depends what is the "home" position of the ellipsoid before rotation).



          Then, the easiest way is to use the invariant properties of the matrix to construct the axis and angle of rotation.



          The trace of the matrix (sum of diagonal) always equals $2costheta+1$, so this is the cheapest way of getting the rotation angle. Note that you still don't know if it's clockwise or anticlockwise (cosine destroys the sign), but you didn't choose the axis orientation either, so it's fine.



          The axis is easiest to get from the antisymmetric part of the matrix. If your matrix is $R$, then calculate



          $$
          R-R^T=
          beginbmatrix
          0&c&-b\ -c & 0 &a \ b & -a & 0
          endbmatrix
          $$
          So, taking the difference of the matrix to its transpose, you can read the axis $(a,b,c)$ from the off-diagonal elements. Essentially, you have
          $$a=R_23-R_32,quad b=R_31-R_13,quad c=R_12-R_21$$



          The length of this vector is not unit, but $2sintheta$, so you must normalize it (you can also verify this angle if it matches the previous one calculated from the trace).



          Now that you have the normalized $(a,b,c)$ axis (unit vector) and the angle $theta$ you can construct the quaternion.



          There's always a possibility that there is another minus that we lost when defining the angle and axis direction - try both and see if it helps.






          share|cite|improve this answer





























            up vote
            0
            down vote













            The phrases you quote seem to be from the LAMMPS manual,
            so it is possible that you will be repeating this calculation for a large number of ellipsoids
            in a molecular dynamics simulation.
            If that's so,
            I feel that it may be helpful to provide an extended version of the excellent answer of
            @orion.
            This is more numerically robust against the danger of dividing by a number close to zero.
            It goes back to S. W. Shepperd,
            "Quaternion from rotation matrix",
            J. Guidance and Control, 1, 223-224 (1978).
            The paper itself is tricky to find online,
            but if you search for "shepperd quaternion" you can turn up several places where the method is described.
            But anyway, here it is.



            There is no need to introduce any trigonometric functions
            (unless you really need to know $theta$, which I doubt).
            Instead of your $(q_w,q_x,q_y,q_z)$ I am going to use the notation $(q_0,q_1,q_2,q_3)$.
            The rotation matrix takes the form
            (which you can find in books such as Classical Mechanics by H Goldstein)
            $$
            beginpmatrix
            R_11 & R_12 & R_13 \
            R_21 & R_22 & R_23 \
            R_31 & R_32 & R_33
            endpmatrix
            =
            beginpmatrix
            q_0^2+q_1^2-q_2^2-q_3^2 & 2(q_1q_2+q_0q_3) & 2(q_1q_3-q_0q_2) \
            2(q_1q_2-q_0q_3) & q_0^2-q_1^2+q_2^2-q_3^2 & 2(q_2q_3+q_0q_1) \
            2(q_1q_3+q_0q_2) & 2(q_2q_3-q_0q_1) & q_0^2-q_1^2-q_2^2+q_3^2
            endpmatrix
            $$
            (It is possible that LAMMPS will have used the transpose of this matrix:
            if so, the following algebra is easily modified).



            Start by calculating the rotation matrix from the unit vectors of the ellipsoid,
            exactly as advised in the answer of @orion.



            Let the trace of the matrix be $T = R_11 + R_22 + R_33$.
            Remember that the quaternions
            are normalized: $q_0^2+q_1^2+q_2^2+q_3^2 = 1$.
            Then, looking at the diagonal elements of $mathbfR$,
            beginalign*
            4q_0^2 & = 1 + 2T - T qquad = 1+T \
            4q_1^2 & = 1 + 2R_11 - T \
            4q_2^2 & = 1 + 2R_22 - T \
            4q_3^2 & = 1 + 2R_33 - T
            endalign*
            The first equation may look a bit odd, but it is written this way
            to make the following choice clear:
            if you choose the largest (most positive) of $ T,R_11,R_22,R_33 $,
            then this clearly corresponds to the largest of $ q_0^2, q_1^2, q_2^2, q_3^2 $.
            Make that choice:
            calculate that component $q_i$ from the corresponding equation.
            You will need to take a square root: you can define this component $q_i$ to be positive
            without loss of generality.



            Now look at the off-diagonal elements of $mathbfR$:
            beginalign*
            4q_0q_1 &= R_23 - R_32 &
            4q_2q_3 &= R_23 + R_32 \
            4q_0q_2 &= R_31 - R_13 &
            4q_3q_1 &= R_31 + R_13 \
            4q_0q_3 &= R_12 - R_21 &
            4q_1q_2 &= R_12 + R_21 .
            endalign*
            Exactly three of these six equations will involve our chosen component $q_i$
            (the one whose magnitude we already know is largest).
            Use these three equations to calculate the other components:
            this will involve dividing by $4q_i$.



            If it turns out that the trace $T$ is larger than any of the diagonal elements of $mathbfR$,
            then $q_0$ will be the largest (in magnitude) of the quaternion parameters,
            and can be calculated as $sqrt(1+T)/4$;
            the components $q_1$, $q_2$, $q_3$ follow
            (from the three equations on the left of
            the above set of six)
            and the procedure reduces to the answer of @orion
            (once the functions involving $theta$ have been eliminated).
            However,
            Shepperd's algorithm guarantees that we avoid numerical troubles if $q_0$ turns out to be small.






            share|cite|improve this answer






















            • Hi @LonelyProf, I did try both yours and orion's methods to calculate the quaternions (and check it several times) but still haven't got a correct match. It is correct that I am doing MD simulations using Lammps. When I use a visualisation tool to displace the atomistic model and the ellipsoid model they don't agree with each other. I calculated the principal axes by taking the eigen-vectors of the moment of inertia tensor, and I'm pretty sure I've got them right because they agree with the molecule coordinates when I plotted them out on a 3D graph... I'm very confused...
              – Nguyen Thi Lan Huong
              Aug 22 at 9:17










            • I'm worried that the comments section may develop into a long question-and-answer thread which is not how maths stackexchange is intended to work. Also, I suspect that the real problem is no longer the one you had at the beginning: I believe that you have two satisfactory answers to that question. It looks more like your problem is computational: either "how to plot an ellipsoid" or "how to understand how LAMMPS works". I am not even sure that such questions would be welcome on scicomp.stackexchange.com because they will not help much with software packages.
              – LonelyProf
              Aug 22 at 9:38










            • Since you have insufficient reputation for chat, let me risk a question or two here. You plotted the inertia axes, along with the original atomistic coordinates, and they look right? But when you try and plot an ellipsoid with those same principal axes, it looks wrong? How are you trying to plot the ellipsoid? Does it rely on knowing the quaternions?
              – LonelyProf
              Aug 22 at 9:46










            • I plotted the atom coordinates and the principal axes using some python code: input coordinates, centre of mass (COM), axis vectors, draw a surface based on atom coordinates, draw lines from COM along the axis vectors, and the axes align right with the surface. For the ellipsoid model, I wrote a Lammps data file with calculated COM and quaternions, ran Lammps for 0 step to output a dump file and visualise it + the atomistic model using OVITO. Sorry about leading this to a long conversation about something else... I also think the problem is how they define ellipsoid orientation in Lammps.
              – Nguyen Thi Lan Huong
              Aug 22 at 10:10










            • I suspect that I can't help much more. If you have tried the transpose of the rotation matrix (or equivalently switching $q_1,q_2,q_3 rightarrow -q_1,-q_2,-q_3$ leaving $q_0$ unchanged), then there's not much else that can be wrong. You can also try plotting the ellipsoid surface by hand using its definition here en.wikipedia.org/wiki/Ellipsoid in the section "As quadric" (defining a matrix $mathbfA$ from your axes: no need for quaternions).
              – LonelyProf
              Aug 22 at 10:29











            Your Answer




            StackExchange.ifUsing("editor", function ()
            return StackExchange.using("mathjaxEditing", function ()
            StackExchange.MarkdownEditor.creationCallbacks.add(function (editor, postfix)
            StackExchange.mathjaxEditing.prepareWmdForMathJax(editor, postfix, [["$", "$"], ["\\(","\\)"]]);
            );
            );
            , "mathjax-editing");

            StackExchange.ready(function()
            var channelOptions =
            tags: "".split(" "),
            id: "69"
            ;
            initTagRenderer("".split(" "), "".split(" "), channelOptions);

            StackExchange.using("externalEditor", function()
            // Have to fire editor after snippets, if snippets enabled
            if (StackExchange.settings.snippets.snippetsEnabled)
            StackExchange.using("snippets", function()
            createEditor();
            );

            else
            createEditor();

            );

            function createEditor()
            StackExchange.prepareEditor(
            heartbeatType: 'answer',
            convertImagesToLinks: true,
            noModals: false,
            showLowRepImageUploadWarning: true,
            reputationToPostImages: 10,
            bindNavPrevention: true,
            postfix: "",
            noCode: true, onDemand: true,
            discardSelector: ".discard-answer"
            ,immediatelyShowMarkdownHelp:true
            );



            );








             

            draft saved


            draft discarded


















            StackExchange.ready(
            function ()
            StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fmath.stackexchange.com%2fquestions%2f2889712%2fhow-to-calculate-quaternions-from-principal-axes-of-an-ellipsoid%23new-answer', 'question_page');

            );

            Post as a guest






























            2 Answers
            2






            active

            oldest

            votes








            2 Answers
            2






            active

            oldest

            votes









            active

            oldest

            votes






            active

            oldest

            votes








            up vote
            1
            down vote



            accepted










            Knowing the orientation vectors is essentially the same as knowing the transformation matrix, which is constructed by putting these vectors (they have to be normalized and orthogonal) into matrix columns.



            First possible issue: arrange the vectors so that they form a positively oriented triad (the determinant of the matrix should be 1, not -1). Otherwise, your matrix is not a rotation but a rotation+reflection. Fixing it is easy, as it's an ellipsoid and you can flip the sign of any main axis vector without changing anything. Second possible issue is that the order of axes is different and you are putting the long axis where the short axis should be (it depends what is the "home" position of the ellipsoid before rotation).



            Then, the easiest way is to use the invariant properties of the matrix to construct the axis and angle of rotation.



            The trace of the matrix (sum of diagonal) always equals $2costheta+1$, so this is the cheapest way of getting the rotation angle. Note that you still don't know if it's clockwise or anticlockwise (cosine destroys the sign), but you didn't choose the axis orientation either, so it's fine.



            The axis is easiest to get from the antisymmetric part of the matrix. If your matrix is $R$, then calculate



            $$
            R-R^T=
            beginbmatrix
            0&c&-b\ -c & 0 &a \ b & -a & 0
            endbmatrix
            $$
            So, taking the difference of the matrix to its transpose, you can read the axis $(a,b,c)$ from the off-diagonal elements. Essentially, you have
            $$a=R_23-R_32,quad b=R_31-R_13,quad c=R_12-R_21$$



            The length of this vector is not unit, but $2sintheta$, so you must normalize it (you can also verify this angle if it matches the previous one calculated from the trace).



            Now that you have the normalized $(a,b,c)$ axis (unit vector) and the angle $theta$ you can construct the quaternion.



            There's always a possibility that there is another minus that we lost when defining the angle and axis direction - try both and see if it helps.






            share|cite|improve this answer


























              up vote
              1
              down vote



              accepted










              Knowing the orientation vectors is essentially the same as knowing the transformation matrix, which is constructed by putting these vectors (they have to be normalized and orthogonal) into matrix columns.



              First possible issue: arrange the vectors so that they form a positively oriented triad (the determinant of the matrix should be 1, not -1). Otherwise, your matrix is not a rotation but a rotation+reflection. Fixing it is easy, as it's an ellipsoid and you can flip the sign of any main axis vector without changing anything. Second possible issue is that the order of axes is different and you are putting the long axis where the short axis should be (it depends what is the "home" position of the ellipsoid before rotation).



              Then, the easiest way is to use the invariant properties of the matrix to construct the axis and angle of rotation.



              The trace of the matrix (sum of diagonal) always equals $2costheta+1$, so this is the cheapest way of getting the rotation angle. Note that you still don't know if it's clockwise or anticlockwise (cosine destroys the sign), but you didn't choose the axis orientation either, so it's fine.



              The axis is easiest to get from the antisymmetric part of the matrix. If your matrix is $R$, then calculate



              $$
              R-R^T=
              beginbmatrix
              0&c&-b\ -c & 0 &a \ b & -a & 0
              endbmatrix
              $$
              So, taking the difference of the matrix to its transpose, you can read the axis $(a,b,c)$ from the off-diagonal elements. Essentially, you have
              $$a=R_23-R_32,quad b=R_31-R_13,quad c=R_12-R_21$$



              The length of this vector is not unit, but $2sintheta$, so you must normalize it (you can also verify this angle if it matches the previous one calculated from the trace).



              Now that you have the normalized $(a,b,c)$ axis (unit vector) and the angle $theta$ you can construct the quaternion.



              There's always a possibility that there is another minus that we lost when defining the angle and axis direction - try both and see if it helps.






              share|cite|improve this answer
























                up vote
                1
                down vote



                accepted







                up vote
                1
                down vote



                accepted






                Knowing the orientation vectors is essentially the same as knowing the transformation matrix, which is constructed by putting these vectors (they have to be normalized and orthogonal) into matrix columns.



                First possible issue: arrange the vectors so that they form a positively oriented triad (the determinant of the matrix should be 1, not -1). Otherwise, your matrix is not a rotation but a rotation+reflection. Fixing it is easy, as it's an ellipsoid and you can flip the sign of any main axis vector without changing anything. Second possible issue is that the order of axes is different and you are putting the long axis where the short axis should be (it depends what is the "home" position of the ellipsoid before rotation).



                Then, the easiest way is to use the invariant properties of the matrix to construct the axis and angle of rotation.



                The trace of the matrix (sum of diagonal) always equals $2costheta+1$, so this is the cheapest way of getting the rotation angle. Note that you still don't know if it's clockwise or anticlockwise (cosine destroys the sign), but you didn't choose the axis orientation either, so it's fine.



                The axis is easiest to get from the antisymmetric part of the matrix. If your matrix is $R$, then calculate



                $$
                R-R^T=
                beginbmatrix
                0&c&-b\ -c & 0 &a \ b & -a & 0
                endbmatrix
                $$
                So, taking the difference of the matrix to its transpose, you can read the axis $(a,b,c)$ from the off-diagonal elements. Essentially, you have
                $$a=R_23-R_32,quad b=R_31-R_13,quad c=R_12-R_21$$



                The length of this vector is not unit, but $2sintheta$, so you must normalize it (you can also verify this angle if it matches the previous one calculated from the trace).



                Now that you have the normalized $(a,b,c)$ axis (unit vector) and the angle $theta$ you can construct the quaternion.



                There's always a possibility that there is another minus that we lost when defining the angle and axis direction - try both and see if it helps.






                share|cite|improve this answer














                Knowing the orientation vectors is essentially the same as knowing the transformation matrix, which is constructed by putting these vectors (they have to be normalized and orthogonal) into matrix columns.



                First possible issue: arrange the vectors so that they form a positively oriented triad (the determinant of the matrix should be 1, not -1). Otherwise, your matrix is not a rotation but a rotation+reflection. Fixing it is easy, as it's an ellipsoid and you can flip the sign of any main axis vector without changing anything. Second possible issue is that the order of axes is different and you are putting the long axis where the short axis should be (it depends what is the "home" position of the ellipsoid before rotation).



                Then, the easiest way is to use the invariant properties of the matrix to construct the axis and angle of rotation.



                The trace of the matrix (sum of diagonal) always equals $2costheta+1$, so this is the cheapest way of getting the rotation angle. Note that you still don't know if it's clockwise or anticlockwise (cosine destroys the sign), but you didn't choose the axis orientation either, so it's fine.



                The axis is easiest to get from the antisymmetric part of the matrix. If your matrix is $R$, then calculate



                $$
                R-R^T=
                beginbmatrix
                0&c&-b\ -c & 0 &a \ b & -a & 0
                endbmatrix
                $$
                So, taking the difference of the matrix to its transpose, you can read the axis $(a,b,c)$ from the off-diagonal elements. Essentially, you have
                $$a=R_23-R_32,quad b=R_31-R_13,quad c=R_12-R_21$$



                The length of this vector is not unit, but $2sintheta$, so you must normalize it (you can also verify this angle if it matches the previous one calculated from the trace).



                Now that you have the normalized $(a,b,c)$ axis (unit vector) and the angle $theta$ you can construct the quaternion.



                There's always a possibility that there is another minus that we lost when defining the angle and axis direction - try both and see if it helps.







                share|cite|improve this answer














                share|cite|improve this answer



                share|cite|improve this answer








                edited Aug 21 at 12:16

























                answered Aug 21 at 12:07









                orion

                12k11832




                12k11832




















                    up vote
                    0
                    down vote













                    The phrases you quote seem to be from the LAMMPS manual,
                    so it is possible that you will be repeating this calculation for a large number of ellipsoids
                    in a molecular dynamics simulation.
                    If that's so,
                    I feel that it may be helpful to provide an extended version of the excellent answer of
                    @orion.
                    This is more numerically robust against the danger of dividing by a number close to zero.
                    It goes back to S. W. Shepperd,
                    "Quaternion from rotation matrix",
                    J. Guidance and Control, 1, 223-224 (1978).
                    The paper itself is tricky to find online,
                    but if you search for "shepperd quaternion" you can turn up several places where the method is described.
                    But anyway, here it is.



                    There is no need to introduce any trigonometric functions
                    (unless you really need to know $theta$, which I doubt).
                    Instead of your $(q_w,q_x,q_y,q_z)$ I am going to use the notation $(q_0,q_1,q_2,q_3)$.
                    The rotation matrix takes the form
                    (which you can find in books such as Classical Mechanics by H Goldstein)
                    $$
                    beginpmatrix
                    R_11 & R_12 & R_13 \
                    R_21 & R_22 & R_23 \
                    R_31 & R_32 & R_33
                    endpmatrix
                    =
                    beginpmatrix
                    q_0^2+q_1^2-q_2^2-q_3^2 & 2(q_1q_2+q_0q_3) & 2(q_1q_3-q_0q_2) \
                    2(q_1q_2-q_0q_3) & q_0^2-q_1^2+q_2^2-q_3^2 & 2(q_2q_3+q_0q_1) \
                    2(q_1q_3+q_0q_2) & 2(q_2q_3-q_0q_1) & q_0^2-q_1^2-q_2^2+q_3^2
                    endpmatrix
                    $$
                    (It is possible that LAMMPS will have used the transpose of this matrix:
                    if so, the following algebra is easily modified).



                    Start by calculating the rotation matrix from the unit vectors of the ellipsoid,
                    exactly as advised in the answer of @orion.



                    Let the trace of the matrix be $T = R_11 + R_22 + R_33$.
                    Remember that the quaternions
                    are normalized: $q_0^2+q_1^2+q_2^2+q_3^2 = 1$.
                    Then, looking at the diagonal elements of $mathbfR$,
                    beginalign*
                    4q_0^2 & = 1 + 2T - T qquad = 1+T \
                    4q_1^2 & = 1 + 2R_11 - T \
                    4q_2^2 & = 1 + 2R_22 - T \
                    4q_3^2 & = 1 + 2R_33 - T
                    endalign*
                    The first equation may look a bit odd, but it is written this way
                    to make the following choice clear:
                    if you choose the largest (most positive) of $ T,R_11,R_22,R_33 $,
                    then this clearly corresponds to the largest of $ q_0^2, q_1^2, q_2^2, q_3^2 $.
                    Make that choice:
                    calculate that component $q_i$ from the corresponding equation.
                    You will need to take a square root: you can define this component $q_i$ to be positive
                    without loss of generality.



                    Now look at the off-diagonal elements of $mathbfR$:
                    beginalign*
                    4q_0q_1 &= R_23 - R_32 &
                    4q_2q_3 &= R_23 + R_32 \
                    4q_0q_2 &= R_31 - R_13 &
                    4q_3q_1 &= R_31 + R_13 \
                    4q_0q_3 &= R_12 - R_21 &
                    4q_1q_2 &= R_12 + R_21 .
                    endalign*
                    Exactly three of these six equations will involve our chosen component $q_i$
                    (the one whose magnitude we already know is largest).
                    Use these three equations to calculate the other components:
                    this will involve dividing by $4q_i$.



                    If it turns out that the trace $T$ is larger than any of the diagonal elements of $mathbfR$,
                    then $q_0$ will be the largest (in magnitude) of the quaternion parameters,
                    and can be calculated as $sqrt(1+T)/4$;
                    the components $q_1$, $q_2$, $q_3$ follow
                    (from the three equations on the left of
                    the above set of six)
                    and the procedure reduces to the answer of @orion
                    (once the functions involving $theta$ have been eliminated).
                    However,
                    Shepperd's algorithm guarantees that we avoid numerical troubles if $q_0$ turns out to be small.






                    share|cite|improve this answer






















                    • Hi @LonelyProf, I did try both yours and orion's methods to calculate the quaternions (and check it several times) but still haven't got a correct match. It is correct that I am doing MD simulations using Lammps. When I use a visualisation tool to displace the atomistic model and the ellipsoid model they don't agree with each other. I calculated the principal axes by taking the eigen-vectors of the moment of inertia tensor, and I'm pretty sure I've got them right because they agree with the molecule coordinates when I plotted them out on a 3D graph... I'm very confused...
                      – Nguyen Thi Lan Huong
                      Aug 22 at 9:17










                    • I'm worried that the comments section may develop into a long question-and-answer thread which is not how maths stackexchange is intended to work. Also, I suspect that the real problem is no longer the one you had at the beginning: I believe that you have two satisfactory answers to that question. It looks more like your problem is computational: either "how to plot an ellipsoid" or "how to understand how LAMMPS works". I am not even sure that such questions would be welcome on scicomp.stackexchange.com because they will not help much with software packages.
                      – LonelyProf
                      Aug 22 at 9:38










                    • Since you have insufficient reputation for chat, let me risk a question or two here. You plotted the inertia axes, along with the original atomistic coordinates, and they look right? But when you try and plot an ellipsoid with those same principal axes, it looks wrong? How are you trying to plot the ellipsoid? Does it rely on knowing the quaternions?
                      – LonelyProf
                      Aug 22 at 9:46










                    • I plotted the atom coordinates and the principal axes using some python code: input coordinates, centre of mass (COM), axis vectors, draw a surface based on atom coordinates, draw lines from COM along the axis vectors, and the axes align right with the surface. For the ellipsoid model, I wrote a Lammps data file with calculated COM and quaternions, ran Lammps for 0 step to output a dump file and visualise it + the atomistic model using OVITO. Sorry about leading this to a long conversation about something else... I also think the problem is how they define ellipsoid orientation in Lammps.
                      – Nguyen Thi Lan Huong
                      Aug 22 at 10:10










                    • I suspect that I can't help much more. If you have tried the transpose of the rotation matrix (or equivalently switching $q_1,q_2,q_3 rightarrow -q_1,-q_2,-q_3$ leaving $q_0$ unchanged), then there's not much else that can be wrong. You can also try plotting the ellipsoid surface by hand using its definition here en.wikipedia.org/wiki/Ellipsoid in the section "As quadric" (defining a matrix $mathbfA$ from your axes: no need for quaternions).
                      – LonelyProf
                      Aug 22 at 10:29















                    up vote
                    0
                    down vote













                    The phrases you quote seem to be from the LAMMPS manual,
                    so it is possible that you will be repeating this calculation for a large number of ellipsoids
                    in a molecular dynamics simulation.
                    If that's so,
                    I feel that it may be helpful to provide an extended version of the excellent answer of
                    @orion.
                    This is more numerically robust against the danger of dividing by a number close to zero.
                    It goes back to S. W. Shepperd,
                    "Quaternion from rotation matrix",
                    J. Guidance and Control, 1, 223-224 (1978).
                    The paper itself is tricky to find online,
                    but if you search for "shepperd quaternion" you can turn up several places where the method is described.
                    But anyway, here it is.



                    There is no need to introduce any trigonometric functions
                    (unless you really need to know $theta$, which I doubt).
                    Instead of your $(q_w,q_x,q_y,q_z)$ I am going to use the notation $(q_0,q_1,q_2,q_3)$.
                    The rotation matrix takes the form
                    (which you can find in books such as Classical Mechanics by H Goldstein)
                    $$
                    beginpmatrix
                    R_11 & R_12 & R_13 \
                    R_21 & R_22 & R_23 \
                    R_31 & R_32 & R_33
                    endpmatrix
                    =
                    beginpmatrix
                    q_0^2+q_1^2-q_2^2-q_3^2 & 2(q_1q_2+q_0q_3) & 2(q_1q_3-q_0q_2) \
                    2(q_1q_2-q_0q_3) & q_0^2-q_1^2+q_2^2-q_3^2 & 2(q_2q_3+q_0q_1) \
                    2(q_1q_3+q_0q_2) & 2(q_2q_3-q_0q_1) & q_0^2-q_1^2-q_2^2+q_3^2
                    endpmatrix
                    $$
                    (It is possible that LAMMPS will have used the transpose of this matrix:
                    if so, the following algebra is easily modified).



                    Start by calculating the rotation matrix from the unit vectors of the ellipsoid,
                    exactly as advised in the answer of @orion.



                    Let the trace of the matrix be $T = R_11 + R_22 + R_33$.
                    Remember that the quaternions
                    are normalized: $q_0^2+q_1^2+q_2^2+q_3^2 = 1$.
                    Then, looking at the diagonal elements of $mathbfR$,
                    beginalign*
                    4q_0^2 & = 1 + 2T - T qquad = 1+T \
                    4q_1^2 & = 1 + 2R_11 - T \
                    4q_2^2 & = 1 + 2R_22 - T \
                    4q_3^2 & = 1 + 2R_33 - T
                    endalign*
                    The first equation may look a bit odd, but it is written this way
                    to make the following choice clear:
                    if you choose the largest (most positive) of $ T,R_11,R_22,R_33 $,
                    then this clearly corresponds to the largest of $ q_0^2, q_1^2, q_2^2, q_3^2 $.
                    Make that choice:
                    calculate that component $q_i$ from the corresponding equation.
                    You will need to take a square root: you can define this component $q_i$ to be positive
                    without loss of generality.



                    Now look at the off-diagonal elements of $mathbfR$:
                    beginalign*
                    4q_0q_1 &= R_23 - R_32 &
                    4q_2q_3 &= R_23 + R_32 \
                    4q_0q_2 &= R_31 - R_13 &
                    4q_3q_1 &= R_31 + R_13 \
                    4q_0q_3 &= R_12 - R_21 &
                    4q_1q_2 &= R_12 + R_21 .
                    endalign*
                    Exactly three of these six equations will involve our chosen component $q_i$
                    (the one whose magnitude we already know is largest).
                    Use these three equations to calculate the other components:
                    this will involve dividing by $4q_i$.



                    If it turns out that the trace $T$ is larger than any of the diagonal elements of $mathbfR$,
                    then $q_0$ will be the largest (in magnitude) of the quaternion parameters,
                    and can be calculated as $sqrt(1+T)/4$;
                    the components $q_1$, $q_2$, $q_3$ follow
                    (from the three equations on the left of
                    the above set of six)
                    and the procedure reduces to the answer of @orion
                    (once the functions involving $theta$ have been eliminated).
                    However,
                    Shepperd's algorithm guarantees that we avoid numerical troubles if $q_0$ turns out to be small.






                    share|cite|improve this answer






















                    • Hi @LonelyProf, I did try both yours and orion's methods to calculate the quaternions (and check it several times) but still haven't got a correct match. It is correct that I am doing MD simulations using Lammps. When I use a visualisation tool to displace the atomistic model and the ellipsoid model they don't agree with each other. I calculated the principal axes by taking the eigen-vectors of the moment of inertia tensor, and I'm pretty sure I've got them right because they agree with the molecule coordinates when I plotted them out on a 3D graph... I'm very confused...
                      – Nguyen Thi Lan Huong
                      Aug 22 at 9:17










                    • I'm worried that the comments section may develop into a long question-and-answer thread which is not how maths stackexchange is intended to work. Also, I suspect that the real problem is no longer the one you had at the beginning: I believe that you have two satisfactory answers to that question. It looks more like your problem is computational: either "how to plot an ellipsoid" or "how to understand how LAMMPS works". I am not even sure that such questions would be welcome on scicomp.stackexchange.com because they will not help much with software packages.
                      – LonelyProf
                      Aug 22 at 9:38










                    • Since you have insufficient reputation for chat, let me risk a question or two here. You plotted the inertia axes, along with the original atomistic coordinates, and they look right? But when you try and plot an ellipsoid with those same principal axes, it looks wrong? How are you trying to plot the ellipsoid? Does it rely on knowing the quaternions?
                      – LonelyProf
                      Aug 22 at 9:46










                    • I plotted the atom coordinates and the principal axes using some python code: input coordinates, centre of mass (COM), axis vectors, draw a surface based on atom coordinates, draw lines from COM along the axis vectors, and the axes align right with the surface. For the ellipsoid model, I wrote a Lammps data file with calculated COM and quaternions, ran Lammps for 0 step to output a dump file and visualise it + the atomistic model using OVITO. Sorry about leading this to a long conversation about something else... I also think the problem is how they define ellipsoid orientation in Lammps.
                      – Nguyen Thi Lan Huong
                      Aug 22 at 10:10










                    • I suspect that I can't help much more. If you have tried the transpose of the rotation matrix (or equivalently switching $q_1,q_2,q_3 rightarrow -q_1,-q_2,-q_3$ leaving $q_0$ unchanged), then there's not much else that can be wrong. You can also try plotting the ellipsoid surface by hand using its definition here en.wikipedia.org/wiki/Ellipsoid in the section "As quadric" (defining a matrix $mathbfA$ from your axes: no need for quaternions).
                      – LonelyProf
                      Aug 22 at 10:29













                    up vote
                    0
                    down vote










                    up vote
                    0
                    down vote









                    The phrases you quote seem to be from the LAMMPS manual,
                    so it is possible that you will be repeating this calculation for a large number of ellipsoids
                    in a molecular dynamics simulation.
                    If that's so,
                    I feel that it may be helpful to provide an extended version of the excellent answer of
                    @orion.
                    This is more numerically robust against the danger of dividing by a number close to zero.
                    It goes back to S. W. Shepperd,
                    "Quaternion from rotation matrix",
                    J. Guidance and Control, 1, 223-224 (1978).
                    The paper itself is tricky to find online,
                    but if you search for "shepperd quaternion" you can turn up several places where the method is described.
                    But anyway, here it is.



                    There is no need to introduce any trigonometric functions
                    (unless you really need to know $theta$, which I doubt).
                    Instead of your $(q_w,q_x,q_y,q_z)$ I am going to use the notation $(q_0,q_1,q_2,q_3)$.
                    The rotation matrix takes the form
                    (which you can find in books such as Classical Mechanics by H Goldstein)
                    $$
                    beginpmatrix
                    R_11 & R_12 & R_13 \
                    R_21 & R_22 & R_23 \
                    R_31 & R_32 & R_33
                    endpmatrix
                    =
                    beginpmatrix
                    q_0^2+q_1^2-q_2^2-q_3^2 & 2(q_1q_2+q_0q_3) & 2(q_1q_3-q_0q_2) \
                    2(q_1q_2-q_0q_3) & q_0^2-q_1^2+q_2^2-q_3^2 & 2(q_2q_3+q_0q_1) \
                    2(q_1q_3+q_0q_2) & 2(q_2q_3-q_0q_1) & q_0^2-q_1^2-q_2^2+q_3^2
                    endpmatrix
                    $$
                    (It is possible that LAMMPS will have used the transpose of this matrix:
                    if so, the following algebra is easily modified).



                    Start by calculating the rotation matrix from the unit vectors of the ellipsoid,
                    exactly as advised in the answer of @orion.



                    Let the trace of the matrix be $T = R_11 + R_22 + R_33$.
                    Remember that the quaternions
                    are normalized: $q_0^2+q_1^2+q_2^2+q_3^2 = 1$.
                    Then, looking at the diagonal elements of $mathbfR$,
                    beginalign*
                    4q_0^2 & = 1 + 2T - T qquad = 1+T \
                    4q_1^2 & = 1 + 2R_11 - T \
                    4q_2^2 & = 1 + 2R_22 - T \
                    4q_3^2 & = 1 + 2R_33 - T
                    endalign*
                    The first equation may look a bit odd, but it is written this way
                    to make the following choice clear:
                    if you choose the largest (most positive) of $ T,R_11,R_22,R_33 $,
                    then this clearly corresponds to the largest of $ q_0^2, q_1^2, q_2^2, q_3^2 $.
                    Make that choice:
                    calculate that component $q_i$ from the corresponding equation.
                    You will need to take a square root: you can define this component $q_i$ to be positive
                    without loss of generality.



                    Now look at the off-diagonal elements of $mathbfR$:
                    beginalign*
                    4q_0q_1 &= R_23 - R_32 &
                    4q_2q_3 &= R_23 + R_32 \
                    4q_0q_2 &= R_31 - R_13 &
                    4q_3q_1 &= R_31 + R_13 \
                    4q_0q_3 &= R_12 - R_21 &
                    4q_1q_2 &= R_12 + R_21 .
                    endalign*
                    Exactly three of these six equations will involve our chosen component $q_i$
                    (the one whose magnitude we already know is largest).
                    Use these three equations to calculate the other components:
                    this will involve dividing by $4q_i$.



                    If it turns out that the trace $T$ is larger than any of the diagonal elements of $mathbfR$,
                    then $q_0$ will be the largest (in magnitude) of the quaternion parameters,
                    and can be calculated as $sqrt(1+T)/4$;
                    the components $q_1$, $q_2$, $q_3$ follow
                    (from the three equations on the left of
                    the above set of six)
                    and the procedure reduces to the answer of @orion
                    (once the functions involving $theta$ have been eliminated).
                    However,
                    Shepperd's algorithm guarantees that we avoid numerical troubles if $q_0$ turns out to be small.






                    share|cite|improve this answer














                    The phrases you quote seem to be from the LAMMPS manual,
                    so it is possible that you will be repeating this calculation for a large number of ellipsoids
                    in a molecular dynamics simulation.
                    If that's so,
                    I feel that it may be helpful to provide an extended version of the excellent answer of
                    @orion.
                    This is more numerically robust against the danger of dividing by a number close to zero.
                    It goes back to S. W. Shepperd,
                    "Quaternion from rotation matrix",
                    J. Guidance and Control, 1, 223-224 (1978).
                    The paper itself is tricky to find online,
                    but if you search for "shepperd quaternion" you can turn up several places where the method is described.
                    But anyway, here it is.



                    There is no need to introduce any trigonometric functions
                    (unless you really need to know $theta$, which I doubt).
                    Instead of your $(q_w,q_x,q_y,q_z)$ I am going to use the notation $(q_0,q_1,q_2,q_3)$.
                    The rotation matrix takes the form
                    (which you can find in books such as Classical Mechanics by H Goldstein)
                    $$
                    beginpmatrix
                    R_11 & R_12 & R_13 \
                    R_21 & R_22 & R_23 \
                    R_31 & R_32 & R_33
                    endpmatrix
                    =
                    beginpmatrix
                    q_0^2+q_1^2-q_2^2-q_3^2 & 2(q_1q_2+q_0q_3) & 2(q_1q_3-q_0q_2) \
                    2(q_1q_2-q_0q_3) & q_0^2-q_1^2+q_2^2-q_3^2 & 2(q_2q_3+q_0q_1) \
                    2(q_1q_3+q_0q_2) & 2(q_2q_3-q_0q_1) & q_0^2-q_1^2-q_2^2+q_3^2
                    endpmatrix
                    $$
                    (It is possible that LAMMPS will have used the transpose of this matrix:
                    if so, the following algebra is easily modified).



                    Start by calculating the rotation matrix from the unit vectors of the ellipsoid,
                    exactly as advised in the answer of @orion.



                    Let the trace of the matrix be $T = R_11 + R_22 + R_33$.
                    Remember that the quaternions
                    are normalized: $q_0^2+q_1^2+q_2^2+q_3^2 = 1$.
                    Then, looking at the diagonal elements of $mathbfR$,
                    beginalign*
                    4q_0^2 & = 1 + 2T - T qquad = 1+T \
                    4q_1^2 & = 1 + 2R_11 - T \
                    4q_2^2 & = 1 + 2R_22 - T \
                    4q_3^2 & = 1 + 2R_33 - T
                    endalign*
                    The first equation may look a bit odd, but it is written this way
                    to make the following choice clear:
                    if you choose the largest (most positive) of $ T,R_11,R_22,R_33 $,
                    then this clearly corresponds to the largest of $ q_0^2, q_1^2, q_2^2, q_3^2 $.
                    Make that choice:
                    calculate that component $q_i$ from the corresponding equation.
                    You will need to take a square root: you can define this component $q_i$ to be positive
                    without loss of generality.



                    Now look at the off-diagonal elements of $mathbfR$:
                    beginalign*
                    4q_0q_1 &= R_23 - R_32 &
                    4q_2q_3 &= R_23 + R_32 \
                    4q_0q_2 &= R_31 - R_13 &
                    4q_3q_1 &= R_31 + R_13 \
                    4q_0q_3 &= R_12 - R_21 &
                    4q_1q_2 &= R_12 + R_21 .
                    endalign*
                    Exactly three of these six equations will involve our chosen component $q_i$
                    (the one whose magnitude we already know is largest).
                    Use these three equations to calculate the other components:
                    this will involve dividing by $4q_i$.



                    If it turns out that the trace $T$ is larger than any of the diagonal elements of $mathbfR$,
                    then $q_0$ will be the largest (in magnitude) of the quaternion parameters,
                    and can be calculated as $sqrt(1+T)/4$;
                    the components $q_1$, $q_2$, $q_3$ follow
                    (from the three equations on the left of
                    the above set of six)
                    and the procedure reduces to the answer of @orion
                    (once the functions involving $theta$ have been eliminated).
                    However,
                    Shepperd's algorithm guarantees that we avoid numerical troubles if $q_0$ turns out to be small.







                    share|cite|improve this answer














                    share|cite|improve this answer



                    share|cite|improve this answer








                    edited Aug 21 at 15:39

























                    answered Aug 21 at 15:29









                    LonelyProf

                    2216




                    2216











                    • Hi @LonelyProf, I did try both yours and orion's methods to calculate the quaternions (and check it several times) but still haven't got a correct match. It is correct that I am doing MD simulations using Lammps. When I use a visualisation tool to displace the atomistic model and the ellipsoid model they don't agree with each other. I calculated the principal axes by taking the eigen-vectors of the moment of inertia tensor, and I'm pretty sure I've got them right because they agree with the molecule coordinates when I plotted them out on a 3D graph... I'm very confused...
                      – Nguyen Thi Lan Huong
                      Aug 22 at 9:17










                    • I'm worried that the comments section may develop into a long question-and-answer thread which is not how maths stackexchange is intended to work. Also, I suspect that the real problem is no longer the one you had at the beginning: I believe that you have two satisfactory answers to that question. It looks more like your problem is computational: either "how to plot an ellipsoid" or "how to understand how LAMMPS works". I am not even sure that such questions would be welcome on scicomp.stackexchange.com because they will not help much with software packages.
                      – LonelyProf
                      Aug 22 at 9:38










                    • Since you have insufficient reputation for chat, let me risk a question or two here. You plotted the inertia axes, along with the original atomistic coordinates, and they look right? But when you try and plot an ellipsoid with those same principal axes, it looks wrong? How are you trying to plot the ellipsoid? Does it rely on knowing the quaternions?
                      – LonelyProf
                      Aug 22 at 9:46










                    • I plotted the atom coordinates and the principal axes using some python code: input coordinates, centre of mass (COM), axis vectors, draw a surface based on atom coordinates, draw lines from COM along the axis vectors, and the axes align right with the surface. For the ellipsoid model, I wrote a Lammps data file with calculated COM and quaternions, ran Lammps for 0 step to output a dump file and visualise it + the atomistic model using OVITO. Sorry about leading this to a long conversation about something else... I also think the problem is how they define ellipsoid orientation in Lammps.
                      – Nguyen Thi Lan Huong
                      Aug 22 at 10:10










                    • I suspect that I can't help much more. If you have tried the transpose of the rotation matrix (or equivalently switching $q_1,q_2,q_3 rightarrow -q_1,-q_2,-q_3$ leaving $q_0$ unchanged), then there's not much else that can be wrong. You can also try plotting the ellipsoid surface by hand using its definition here en.wikipedia.org/wiki/Ellipsoid in the section "As quadric" (defining a matrix $mathbfA$ from your axes: no need for quaternions).
                      – LonelyProf
                      Aug 22 at 10:29

















                    • Hi @LonelyProf, I did try both yours and orion's methods to calculate the quaternions (and check it several times) but still haven't got a correct match. It is correct that I am doing MD simulations using Lammps. When I use a visualisation tool to displace the atomistic model and the ellipsoid model they don't agree with each other. I calculated the principal axes by taking the eigen-vectors of the moment of inertia tensor, and I'm pretty sure I've got them right because they agree with the molecule coordinates when I plotted them out on a 3D graph... I'm very confused...
                      – Nguyen Thi Lan Huong
                      Aug 22 at 9:17










                    • I'm worried that the comments section may develop into a long question-and-answer thread which is not how maths stackexchange is intended to work. Also, I suspect that the real problem is no longer the one you had at the beginning: I believe that you have two satisfactory answers to that question. It looks more like your problem is computational: either "how to plot an ellipsoid" or "how to understand how LAMMPS works". I am not even sure that such questions would be welcome on scicomp.stackexchange.com because they will not help much with software packages.
                      – LonelyProf
                      Aug 22 at 9:38










                    • Since you have insufficient reputation for chat, let me risk a question or two here. You plotted the inertia axes, along with the original atomistic coordinates, and they look right? But when you try and plot an ellipsoid with those same principal axes, it looks wrong? How are you trying to plot the ellipsoid? Does it rely on knowing the quaternions?
                      – LonelyProf
                      Aug 22 at 9:46










                    • I plotted the atom coordinates and the principal axes using some python code: input coordinates, centre of mass (COM), axis vectors, draw a surface based on atom coordinates, draw lines from COM along the axis vectors, and the axes align right with the surface. For the ellipsoid model, I wrote a Lammps data file with calculated COM and quaternions, ran Lammps for 0 step to output a dump file and visualise it + the atomistic model using OVITO. Sorry about leading this to a long conversation about something else... I also think the problem is how they define ellipsoid orientation in Lammps.
                      – Nguyen Thi Lan Huong
                      Aug 22 at 10:10










                    • I suspect that I can't help much more. If you have tried the transpose of the rotation matrix (or equivalently switching $q_1,q_2,q_3 rightarrow -q_1,-q_2,-q_3$ leaving $q_0$ unchanged), then there's not much else that can be wrong. You can also try plotting the ellipsoid surface by hand using its definition here en.wikipedia.org/wiki/Ellipsoid in the section "As quadric" (defining a matrix $mathbfA$ from your axes: no need for quaternions).
                      – LonelyProf
                      Aug 22 at 10:29
















                    Hi @LonelyProf, I did try both yours and orion's methods to calculate the quaternions (and check it several times) but still haven't got a correct match. It is correct that I am doing MD simulations using Lammps. When I use a visualisation tool to displace the atomistic model and the ellipsoid model they don't agree with each other. I calculated the principal axes by taking the eigen-vectors of the moment of inertia tensor, and I'm pretty sure I've got them right because they agree with the molecule coordinates when I plotted them out on a 3D graph... I'm very confused...
                    – Nguyen Thi Lan Huong
                    Aug 22 at 9:17




                    Hi @LonelyProf, I did try both yours and orion's methods to calculate the quaternions (and check it several times) but still haven't got a correct match. It is correct that I am doing MD simulations using Lammps. When I use a visualisation tool to displace the atomistic model and the ellipsoid model they don't agree with each other. I calculated the principal axes by taking the eigen-vectors of the moment of inertia tensor, and I'm pretty sure I've got them right because they agree with the molecule coordinates when I plotted them out on a 3D graph... I'm very confused...
                    – Nguyen Thi Lan Huong
                    Aug 22 at 9:17












                    I'm worried that the comments section may develop into a long question-and-answer thread which is not how maths stackexchange is intended to work. Also, I suspect that the real problem is no longer the one you had at the beginning: I believe that you have two satisfactory answers to that question. It looks more like your problem is computational: either "how to plot an ellipsoid" or "how to understand how LAMMPS works". I am not even sure that such questions would be welcome on scicomp.stackexchange.com because they will not help much with software packages.
                    – LonelyProf
                    Aug 22 at 9:38




                    I'm worried that the comments section may develop into a long question-and-answer thread which is not how maths stackexchange is intended to work. Also, I suspect that the real problem is no longer the one you had at the beginning: I believe that you have two satisfactory answers to that question. It looks more like your problem is computational: either "how to plot an ellipsoid" or "how to understand how LAMMPS works". I am not even sure that such questions would be welcome on scicomp.stackexchange.com because they will not help much with software packages.
                    – LonelyProf
                    Aug 22 at 9:38












                    Since you have insufficient reputation for chat, let me risk a question or two here. You plotted the inertia axes, along with the original atomistic coordinates, and they look right? But when you try and plot an ellipsoid with those same principal axes, it looks wrong? How are you trying to plot the ellipsoid? Does it rely on knowing the quaternions?
                    – LonelyProf
                    Aug 22 at 9:46




                    Since you have insufficient reputation for chat, let me risk a question or two here. You plotted the inertia axes, along with the original atomistic coordinates, and they look right? But when you try and plot an ellipsoid with those same principal axes, it looks wrong? How are you trying to plot the ellipsoid? Does it rely on knowing the quaternions?
                    – LonelyProf
                    Aug 22 at 9:46












                    I plotted the atom coordinates and the principal axes using some python code: input coordinates, centre of mass (COM), axis vectors, draw a surface based on atom coordinates, draw lines from COM along the axis vectors, and the axes align right with the surface. For the ellipsoid model, I wrote a Lammps data file with calculated COM and quaternions, ran Lammps for 0 step to output a dump file and visualise it + the atomistic model using OVITO. Sorry about leading this to a long conversation about something else... I also think the problem is how they define ellipsoid orientation in Lammps.
                    – Nguyen Thi Lan Huong
                    Aug 22 at 10:10




                    I plotted the atom coordinates and the principal axes using some python code: input coordinates, centre of mass (COM), axis vectors, draw a surface based on atom coordinates, draw lines from COM along the axis vectors, and the axes align right with the surface. For the ellipsoid model, I wrote a Lammps data file with calculated COM and quaternions, ran Lammps for 0 step to output a dump file and visualise it + the atomistic model using OVITO. Sorry about leading this to a long conversation about something else... I also think the problem is how they define ellipsoid orientation in Lammps.
                    – Nguyen Thi Lan Huong
                    Aug 22 at 10:10












                    I suspect that I can't help much more. If you have tried the transpose of the rotation matrix (or equivalently switching $q_1,q_2,q_3 rightarrow -q_1,-q_2,-q_3$ leaving $q_0$ unchanged), then there's not much else that can be wrong. You can also try plotting the ellipsoid surface by hand using its definition here en.wikipedia.org/wiki/Ellipsoid in the section "As quadric" (defining a matrix $mathbfA$ from your axes: no need for quaternions).
                    – LonelyProf
                    Aug 22 at 10:29





                    I suspect that I can't help much more. If you have tried the transpose of the rotation matrix (or equivalently switching $q_1,q_2,q_3 rightarrow -q_1,-q_2,-q_3$ leaving $q_0$ unchanged), then there's not much else that can be wrong. You can also try plotting the ellipsoid surface by hand using its definition here en.wikipedia.org/wiki/Ellipsoid in the section "As quadric" (defining a matrix $mathbfA$ from your axes: no need for quaternions).
                    – LonelyProf
                    Aug 22 at 10:29













                     

                    draft saved


                    draft discarded


























                     


                    draft saved


                    draft discarded














                    StackExchange.ready(
                    function ()
                    StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fmath.stackexchange.com%2fquestions%2f2889712%2fhow-to-calculate-quaternions-from-principal-axes-of-an-ellipsoid%23new-answer', 'question_page');

                    );

                    Post as a guest













































































                    這個網誌中的熱門文章

                    How to combine Bézier curves to a surface?

                    Mutual Information Always Non-negative

                    Why am i infinitely getting the same tweet with the Twitter Search API?