Some alignment of quaternion formalism with the referenced works

[bogdan-tanygin]

In regards to all the key quaternion related equations in the rotation.cpp... There are references sonnenschein1985 and Omelyan1998 though the code and Martys1999 in the code header and in the documentation. A comparison of the formulae provides a match between these 3 works. For the source code: if one compares each component index of the quaternion and such indices of the torque, omega and other Cartesian (hereinafter, it means 3D space) entities then there is a mismatch with the reference work formulae. Of course, all the formalism works fine now. A quick analysis result: this is just a formal redefinition of the quaternion space with some other alignments: the 4-dim space axes permutation and some inversion. Note, that the Cartesian components are the same. This is why most of formulae looks different. Probably, the real original reference work was different or some other formalism has been used initially. Just one example of practical confusion of this mismatch (there are others as well: with the director, etc):

• from Es. documentation: For a value of (1,0,0,0), the body and space frames coincide.

• from Eq. 3 of sonnenschein1985: see the image (the rotation matrix):
  
  Obviously, the ksi is 0-th element of the quaternion. If we set 1 there (from q=1 0 0 0 value) then we have the rotational matrix diagonal elements: -1 1 -1. The Eq. 3 will require q=0 0 0 1 in this case. Actually, we can see the the Eq. 3 is different compare to the matrix in define_rotation_matrix().

The same for other equations. Nothing serious, but it, probably, leads to some confusions when someone will try to use the quaternions from Es. in a context of other equations and the formalism and simulation methods which are based on the conventions of the Martys1999, other references with such formalism.

Unless, I seriously misunderstood something, please correct me if this is the case... I'm ready to help: either documentation based or the source code adaptations (permute the indices in the Q formulae, maybe align some tests, etc). You can decide, but I suggest the latter solution cause it is really sophisticated to work with the existing quaternion space when trying to develop some tailored extra features of scientific computing... At least for me :-)

[bogdan-tanygin]

a bit related discussions: #1048, #1854

[bogdan-tanygin]

this can be related to the notation convention options: either 1st or last component of the quaternion is real. However, probably, this does not change my above comments ;-)

[RudolfWeeber]

It is not clear to me, whether this would be a user-visible change. Does this change the numeric values of ParticleHandle.quat representing a specific orientation?

[bogdan-tanygin]

I think that it will swap 1st and last quaternion component visible by user. Maybe some signs. I've already seen 2 options of these formalism in different books. Maybe I will just add a reference to the textbook additionally to the present references with a note of this difference? Sincerely yours / mit freundlichen Grüßen, Bogdan Tanygin

…

On Thu, 12 Sep 2019, 12:32 RudolfWeeber, ***@***.***> wrote: It is not clear to me, whether this would be a user-visible change. Does this change the numeric values of ParticleHandle.quat representing a specific orientation? — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <#2964?email_source=notifications&email_token=AC43HUK7SBS7WPR77NQPJWTQJIEEPA5CNFSM4H4VVON2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOD6RIWFI#issuecomment-530746133>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AC43HUK7FR4ETP65FIJX653QJIEEPANCNFSM4H4VVONQ> .

[pdebuyl]

Hey, just pitching in :-)

One important issue is whether one can store quaternions in the output data. If so, I would like to support -- if possible -- to conform to the storage scheme of a nice quaternion library in Python, rowan. Their scheme is "scalar i j k" components.

For processing the rotational dynamics via, say, H5MD files, this would make it very convenient :-)

[bogdan-tanygin]

in this case, I suppose, that existing quaternion scheme is already fine. Cause, for a value of (1,0,0,0), the body and space frames coincide. I need to check, but I think so. So, only a reference correction (addition), perhaps, is needed. I can do it if needed.

[jngrad]

In my opinion, the current formalism is espresso is not satisfactory. The indices do not match those from other libraries like boost::math::quaternion, and operations like rotation inversion and first/second derivative would probably be easier to maintain if we had a quaternion class with its own set of quaternion-specific operators.

[RudolfWeeber]

When I implemented the code rotating a single particle, I had the imperssion that Espresso uses quaternions the way, they are described in the Wikipedia.

I have reservations about changing the meaning of quaternion indices in a way that is visible to the user. That would be a really nasty silent interface change. Users specify quaternions directly, e.g., for virtual sites which have a fixed orientation with respect to the orientation of the central particle.

I see the benefit of replacing stuff by library code. On the other hand, I'm not sure the benefits in this particular case warrant the potential hazard a definition change can cause to users.
Particularly, since such a change doesn't provide significant new functionality to the user and the current code is well tested.

[bogdan-tanygin]

In my opinion, the current formalism is espresso is not satisfactory. The indices do not match those from other libraries like boost::math::quaternion, and operations like rotation inversion and first/second derivative would probably be easier to maintain if we had a quaternion class with its own set of quaternion-specific operators.

@jngrad could you please specify which indices do not match? Is the boost quaternion "scalar i j k"?

[bogdan-tanygin]

please correct me if I'm wrong. I think, I've just found answer:
https://cs.brown.edu/~jwicks/boost/libs/math/quaternion/quaternion.html
... real and R_component_1 return the same value. ...
Hence, yes, should be "scalar i j k".

[jngrad]

Nevermind, quaternion<T> & boost::math::quaternion::operator *=(quaternion<X> const & rhs) has the exact same formula as in our multiplication function if I re-arrange the terms properly. Both boost and espresso use the {scalar, vector} aka (w, i, j, k) format.

[bogdan-tanygin]

great! So, after the work on #3164 I suggest to add some reference which align referred works Martys1999/Omelyan1998/sonnenschein1985 towards our implemented formalism. It can be e.g. Allen, Michael P., and Dominic J. Tildesley. Computer simulation of liquids. Oxford university press, 2017. considering highly cited textbooks. I can help with it, probably, this is just comments in the source code + documentation though..