[Feature Proposal] Fix momentum angular for non-rigid groups of atoms

Dear LAMMPS community,

I’m Arnau Jurado Romero and I’m currently working on my PhD at the Politechnic University of Catalonia (UPC). I want to implement and add a new feature to LAMMPS and I need some guidance as to how to organize its components.

As you may know, LAMMPS currently has a fix momentum for a group of atoms which uses the angular momentum and inertia tensor computation methods implemented in the group.cpp code (when called with the ‘angular’ keyword). These methods assume that the group of atoms is rigid in time, i.e. their relative positions don’t change. This results in incorrect angular velocities if the relative positions do change which makes it so the angular momentum of the group is not removed, instead some rotations is left and some vibration is removed (or added). For most cases this is not a problem but I have found that it does cause the group to rotate and drift over long times.

I have done substantial work in codes for analyzing the vibrations of molecules and to accurately describe those vibrations you need to compute the precisely the rotational velocity taking into account relative positions of the atoms of the molecule. More precisely, the Eckart rotation matrix has to be computed, which relates the original equilibrium reference frame of the molecule with the instantaneous reference frame (the Eckart reference frame). With the equilibrium coordinates in this instantaneous reference frame one can compute an angular momentum and inertia tensor that does not assume rigidity.

I have implemented and tested this methodology in FORTRAN and recently I have implemented it in LAMMPS as a fix which is very similar in structure to fix momentum but with the added computation of the Eckart frame, etc. But I have some trouble in how to organize the code to be submitted in github for LAMMPS.

Here I summarize the method so you can see the general requirements:

  1. Store the equilibrium coordinates of the molecule (group of atoms). I have used a fix store/peratom for this.
  2. At the end step, compute a 4x4 matrix and diagonalize it.
  3. Using the eigenvector of the smallest eigenvalue of this 4x4 matrix compute the Eckart rotation matrix.
  4. Rotate the equilibrium coordinates using this rotation matrix.
  5. Use the adapted inertia tensor and angular momentum formulas.
  6. Compute omega and remove angular momentum.

Point 1 makes it so the group this fix is applied to cannot be extended or reduced (in contrast to the regular fix momentum). In my current implementation I just take the coordinates when the fix is defined as ‘equilibrium’ but it would probably be best to be inputted by file.

Point 2 has required to add to math_eigen some code for 4x4 matrices, since the current implementation only supports 3x3.

Point 6 is exactly the same as in the regular fix momentum and for point 5 I have just substituted the calls to group->angmom and group->inertia with my own functions. The nice thing about the adapted formulas is that they reduced to the regular formulas for rigid molecules/groups (when the rotated equilibrium coordinates match with the instantaneous coordinates).

As you can see, you can use this for removing angular momentum but it gives lots of ‘bonus’ results that LAMMPS does not have currently, mainly an accurate angular momentum, inertia moment and angular velocity for non-rigid groups of atoms, which would correspond currently with the variables angmom(), inertia() and omega().

So, my questions are the following:

  • Should I add my modified angmom and inertia methods to group.cpp? If so, how should I deal with the other requirements (storage of eq. coordinates, for example). This seems like a bad idea, since it is a core file in LAMMPS.
  • Should I only implement the modified fix momentum and add a compute to get the other quantities? Or implement only the modified fix momentum and see if some users want access to the other quantities. (I personally don’t have a need for them, but I do have a need for the fix momentum).
  • Can I edit math_eigen.cpp and .h directly and include it in my pull request or do I need to add a different file?

Any suggestions on how to organize the code would be greatly appreciated. As I say, the core of the code for making what I described has already been done and tested by me. My doubts are about how to package these tools for inclusion in LAMMPS.

Of course, when I submit the code to github I will include all the appropriate articles and references. Or you can ask for them, I just don’t want to make this post much longer.

I hope this channel is the appropriate one for this kind of discussion.

You have to provide some proof to back this statement up. From looking at the source code this is not the case. All group properties are recomputed every step the fix is executed.

Could it be that you are misunderstanding what this fix is supposed to be doing?

Of course. I find that section II.A. of J. Chem. Phys. 107 , 1394 (1997) is a good starting point for the separation of rotations and vibrations in molecular motion. There are numerous articles that I can provide (and are cited by this one) that explore more in depth this topic.

In summary, when computing omega with the standard inertia tensor (and angular momentum) some part of this omega is actually vibrational motion. If you try to remove the rotational velocities derived from this omega (omega x r_cm) the vibrational motion will be modified (and you may not even remove all rotation!). To avoid this what should be computed is the angular velocity of the equilibrium coordinates in the reference frame that rotates with the molecule.

Finding this frame is a whole other issue, but it is widely accepted that the frame that satisfies the Eckart conditions (a.k.a the Eckart frame) does the best job at separating vibration and rotation (It minimizes the last term of equation 12 of the article, also known as the Coriolis term, which can be understood as the vibrational-rotational coupling).

So, computing this omega free of vibration requires a slightly different formula: equation 18 and 19 of the article. As you will be able to tell, these formulas reduce to the usual definitions of inertia tensor and angular momentum if the equilibrium coordinates are the same as the actual CoM coordiantes. This is what I mean with the molecule being rigid, it doesn’t vibrate.

Could be. I understand that this fix should stop a group of atoms from rotating or, in other words, make it so their collective rotation energy is 0. For groups of atoms that vibrate (such as non-rigid molecules) the current implementation is not exactly correct.

I realize that this is a very niche issue and that fix momentum was probably not intended to be used like I intend to, that’s why I wanted to implement this.

If you wish for more references and/or further explanations I am happy to provide them.

fix momentum knows nothing about molecules or similar. Its purpose is rather simple: remove the current momentum. Initially, this was only done for translation, rotation support was added later.
In some way it is considered a “hack” because in a system without cutoffs, external forces, and perfect precision, there should be no residual momentum formed if there is none initially.
It is not expected that it will stop rotation or translation entirely only that it avoids a situation, where the residual motion of the fix group will be significant enough to impact analysis or thermostatting. That is why it is crucial to figure out the optimal interval to have not too much of a change made when the fix becomes active but it should not manipulate the system too often so that the (regular) manipulation will trigger some unphysical effects/dynamics, and that first all other options to reduce a drift should be applied first.

If I wanted to (smoothly) stop a system from translating (or rotating) I would prefer applying a sufficiently small restraint force.

I agree that using this fix to remove rotation regularly is unphysical and can cause issues. Still, when applied to some reduced group of atoms (i.e. not the whole system but, for example, a molecule) this fix does not remove only momentum, it also removes vibrations. However, it is clear to me now that this was never the intention of this fix, thank you for that clarification.

I’m beginning to think that maybe presenting this as modification to fix momentum was not a good idea. I still think there is some value to adding a compute for the equilibrium coordinates of the molecule in the reference frame that rotates with the molecule, or maybe the omega of this rotating frame.

I don’t know if there is even demand for this kind of result. I do, at least. Like I said most of the code is already done, it’s just a question of how to package it. What do you think, Axel?

It depends a lot on what the intended application scenario is. I think grouping by molecule (ID) is rather arbitrary. LAMMPS has the “chunk” infrastructure to do flexible grouping (e.g. for clusters, fragments, bins, molecules, etc.), so implementing a fix momentum/chunk may be a more suitable target. We have been trying to avoid adding features of this kind that use other ways of grouping atoms. Changing that would require some additional programming, though.

LAMMPS has a lot of components that are only useful for niche applications, so there is no principal problem to add code for another such applications. This one would go to the EXTRA-FIX package. The most important component would be a good documentation that clearly explains what the use case is and best also one simple example input to demonstrate its use.

The second challenge is to write sufficiently clean and maintainable code. There are some suggestions and requirements in the LAMMPS manual at 3.3. LAMMPS programming style and requirements for contributions — LAMMPS documentation and 4.8. Notes for developers and code maintainers — LAMMPS documentation . Also we are trying to increase internal consistency and utilization of modernized internal APIs. Some information on that is posted here: LAMMPS Code Clinic 2022

The discussion on the details of the implementation is best done on GitHub directly. You can elect to submit a draft pull request and invite comments and suggestions while you are still working on it. Sometimes, showing your code early can get you some feedback that would save you some effort compared to having to make changes later after completing an implementation.

There is already a fix momentum/chunk in EXTRA-FIX, so I will think of some other name. Besides, when I used these tools it was only for single molecule in water so using the standard grouping has been enough for me. I can see the use in implementing a chunk based grouping (if you had more than one of the same molecule in the system) but I will stick to standard grouping for the initial implementation.

Ok, I will prepare a draft to submit to github and see what other developers and users have to say about it. Thank you for the resources and for the help.

Hi Arnau –

Out of curiosity, is this meant to be a per-molecule transformation? That is, is the fix meant to perform its adjustments for each molecule in its own Eckart frame?

I am part of a team independently working on other per-molecule quantities – let me know if you would like to compare notes.


Yes, each molecule has it’s own rotating (Eckart) reference frame so the idea is to apply the transformation on a per-molecule basis. In my research I have only applied this transformation to a single molecule in the system (with a post-processing fortran code that I developed). So that is why I thought of this potential fix/compute for a group of a few atoms. But this could easily be expanded to compute all the Eckart frames of all molecules in the system by using the chunk based grouping, like Axel suggested.

Sure, that would be wonderful. What quantities is your team working on?