Question on the calculation of the Steinhardt Order Parameter

LAMMPS LAMMPS Development
1

Hi,

I am trying to calculate some quantities from MD simulations that require the calculation of spherical harmonics as a function of the distance between two atoms. I am trying to write this as a compute function in LAMMPS. As I do not have much experience in C++ coding, I am following the compute orientorder/atom function in LAMMPS. As I am reading through the orientorder/atom compute (https://github.com/lammps/lammps/blob/35214120ad350b7a7fb6e929f6abf467e21fcf8d/src/compute_orientorder_atom.cpp), I am confused about the use of swap3 function (https://github.com/lammps/lammps/blob/35214120ad350b7a7fb6e929f6abf467e21fcf8d/src/compute_orientorder_atom.cpp#L344).

It seems to me it is a sorting function of some kind. My question is why is there a need to sort the arrays defined in the compute orientorder/atom function before calculating the Q_l and W_l variables (via calc_boop() function) needed for the order parameters? Does it help with performance?

Max

2

This link points to a member function called select3(), not swap3(). There is a macro defined before it called SWAP3() (and also not swap3()).

You can see the purpose of the select3() function from reading the comments where the function is called and where it is defined:

      // if nnn > 0, use only nearest nnn neighbors

      if (nnn > 0) {
        select3(nnn, ncount, distsq, nearest, rlist);
        ncount = nnn;
      }

/* ----------------------------------------------------------------------
   select3 routine from Numerical Recipes (slightly modified)
   find k smallest values in array of length n
   sort auxiliary arrays at same time
------------------------------------------------------------------------- */

[...]
void ComputeOrientOrderAtom::select3(int k, int n, double *arr, int *iarr, double **arr3)
{

The SWAP3() macro just swaps the values of two 3-vectors.

It is a partial sort.

The compute by default only considers the 12 nearest neighbors (or a different number if the “nnn” keyword is used) and thus the calc_boop() computation only needs to be done for those “nnn” values.

3

Hi there,
I am not familiar with your type of problem (so I might be wrong), but I noticed someone already made a code for calculating Steinhardt’s bond orientational order parameters. It might be helpful for you to check it:
https://docs.pyscal.org/en/latest/methods/02_steinhardt.html