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:
- Store the equilibrium coordinates of the molecule (group of atoms). I have used a fix store/peratom for this.
- At the end step, compute a 4x4 matrix and diagonalize it.
- Using the eigenvector of the smallest eigenvalue of this 4x4 matrix compute the Eckart rotation matrix.
- Rotate the equilibrium coordinates using this rotation matrix.
- Use the adapted inertia tensor and angular momentum formulas.
- 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.