Just a comment regarding lammps documentation for compute gyration (and gyration/molecule). Looking at the source code and verified by manual computations from trajectory dumps, it seems that when the compute calculates the actual Rg (for a group or per molecule) it stores the square root of the rg in units of distance. When the computes calculate the gyration tensor it stores the elements of the tensor squared. In the output section of the documentation it describes all output as being in units of length (not length squared for the tensor). It might be worth clarifying the documentation in this regard. At least for me, every time I go back to this command I have to re-read the docs, which never helps. I always have to do mini-tests to make sure whether I am getting Rg or Rg^2. Maybe there is something I am missing?
Perhaps the following changes are enough:
“The 6 components of the vector are ordered xx, yy, zz, xy, xz, yz (each squared)”
“This compute calculates a global scalar (Rg) and a global vector of length 6 (Rg squared tensor)”
“The scalar and vector values will be in length and length^2 units, respectively”
Similar changes for compute/gyration molecule would be needed.
ok - we could do what you suggest, or we could
take the sqrt() of the vector components. Would
the latter be simpler and more consistent? I’m not
sure what the standard tensor formulation should be.
For the tensors, the user might be interested in diagonalizing the matrix to get the principal moments and subsequently get shape descriptors like asphericity, acylindricity, and the relative shape anisotropy. It is my understanding that all these quantities are based on the squares of the tensor elements so in this regard the implementation is fine but inconsistent with the Rg scalar quantity. If consistency is the goal, I believe we can take the sqrt() of the tensor and let the user square the results post-simulation.
done - I made them all sqrt() values and adjusted
the doc pages to make this clearer.
will be in the next path
I thought about this a little more and I have a question. If one wants to use fix ave/time on this tensor with your new patched version, how could one average the squared values? Would we have to define an equal style variable of the tensor squared and use that in the fix ave/time command?
The reason I ask is because people will typically want the <Rxy^2> for example and not the ^2 value.
I think your solution still makes everything more consistent. I should be careful what I wish for, since this will complicate my scripts a little bit and I will have to square the values manually before averaging them.
Randomly parsing stuff here.
It makes more sense to me (to avoid the most confusion) to use the radius of gyration squared and its tensor form and just note it up front . I say this mostly due to its use in calculating rotational kinetic energy, which is why I am guessing anyone would be interested in it(if we’re talking about single molecules/aggregates). Its not much work to sqrt() something. AND there are better structure metrics if you’re talking about things which are not a single molecule/aggregate.
Just my two cents.
Should hold for fractal dimension as well.
Last piece and I’ll let it die, all of the other metrics (e.g. asphericity and acylindricity) are cast in squares of the diagonal elements of the ‘radius of gyration tensor’ as well.
Eric is right. My concern is that it will be confusing for users to use compute GYRATION and get the rg squared. Then again, if you cannot read the docs it is your own problem.