# Bond-order parameter $g_6(r)$

I would like to know the explicit steps to compute the bond correlation function $$g_6(r) = \langle\Psi_6^*(0) \cdot \Psi_6(r)\rangle$$ for colloids in my experiment (2D). I have the position of all my particles, I know how to compute the individual hexagonal order parameter for particle $$k$$ $$\Psi_{6,k} = \sum_{j=1}^N\frac{\exp[6i\theta_{kj}]}{n_j}$$ where $$n_j$$ is the number of neighbors of the particle $$k$$ (usually computed via Voronoi Tessellation that define the neighbor without cutoff) and $$\theta_{kj}$$ is the angle between $$\vec{r}=\vec{r}_j-\vec{r}_k$$ (the bond) and an arbitrary axis (x-axis for example)

But I cannot see how $$\Psi_{6,k}$$ depends on the position to apply the formula for $$g_6(r)$$. How can I compute $$\Psi_{6k}(r)$$? Is it some kind of average or density function?

• Please define $\theta_{kj}$ and $n_j$, and what is the subscript $k$ as well? Dec 2 '21 at 17:29
• physics.stackexchange.com/questions/461153/… may be of use Dec 2 '21 at 17:32
• Rushing out but looking at the link of the physics forum it seems (in the Questions language) k is the index of a particle, and thus you get position dependence via that index Dec 2 '21 at 17:35
• @IanBush I have updated the question. If $\psi_{6k}(r)$ is the order parameter for the particle $k$, I know the value of $\psi_{6k}(r_k)$ but not the value for all $r$. Then, how can I compute $\psi_{6k}(r)$? Dec 2 '21 at 17:41
• See this which might help (or references therein). Looks to me like they compute the value by summing over angles ($\theta$). Dec 2 '21 at 18:10

Assume we are given a configuration of atoms $$\{{\bf r_i}\}$$ and any other data required to characterise the atomic layout, e.g. any lattice vectors if we have periodic boundary conditions in any of the 3 possible directions. With each of the atoms let us associate a property that can be characterised by a scalar quantity $$\{p_i\}$$. For example $$p_i=1 \space \forall i$$ is one possibility, this will lead to something closely related to the radial distribution function. Another possibility is $$p_i=q_i$$ where $$q_i$$ is the charge on atom $$i$$. Another possibility is $$\Psi_{6,k}$$ as defined in the question. It doesn't really matter what, as long as you can calculate a scalar quantity associated with each of the atoms of interest the recipe below will work [1]

So the first step in the process is for the given configuration calculate $$\{p_i\}$$ for the property of interest. Now we want to calculate the average of $$p_i^\dagger . p_j$$ where atom $$i$$ and atom $$j$$ are separated by a given distance $$d$$. In practice on a computer we divide the distance range into a number of different bins, each of which cover the range $$d$$ to $$d + \delta d$$. So the method will be something like the following - but note I have had to make assumptions about how you wish to normalise the accumulated values; you may have to modify that according to what the relevant literature says. Anyway the essence of the method

1. Set to zero two sets of bins, one to accumulate the correlation terms to, one to count how many atom pairs are in each distance range
2. Loop over all the first particles $$i$$
3. Loop over all the second particles $$j$$[2]
4. Calculate $$v=p_i^\dagger . p_j$$
5. Calculate the distance between atom $$i$$ and atom $$j$$
6. Work out which bin that distance corresponds to
7. Add $$v$$ to the appropriate correlation accumulator bin
8. Increment the appropriate atom pair counter bin
9. Finish the loop over $$j$$
10. Finish the loop over $$i$$
11. The average correlation at each distance is the value of the correlation accumulator bin divided by the atom pair counter bin corresponding to the distance of interest

This generates the contribution for a single configuration - the final step is if appropriate to average over a large number of configurations.

[1] For the more mathematically minded I would normally convolute $$\{p_i\}$$ with delta functions at the atoms positions to represent the distribution of the property of interest. This complication is not needed here, but makes the argument more easily generalised to non-point particles

[2] Looping over all atoms in the $$j$$ loop is not the most efficient way to do this, and will be very slow for large systems. If you care about time to solution you should set up a neighbour list for each atom before the loop, and then in the $$j$$ loop only go over those atoms within the maximum distance of interest from the reference atom $$i$$. this reduces the scaling of the algorithm from $$O(N^2)$$ to $$O(N)$$ where $$N$$ is the number of atoms in the system

• I see! Thank you very much for the explanation, now it is more clear for me. Respect to the normalization I think that I just need to divide for the total number of particles to find the correct function. I'm gonna add at the end of my question the code (write in python) that I'm doing for people that maybe need it in the future! Dec 3 '21 at 16:38
• I almost wrote some code, but as a Fortran or C guy I didn't want to frighten people ... Dec 3 '21 at 17:26
• It looks like $p_i$ and $p_j$ become vectors, rather than scalars in step 4. Dec 3 '21 at 17:34
• A dot can be used to mean multiply, it's not just for dot products. I'll remove if it causes confusion Dec 3 '21 at 17:45
• It's fine for me, a presume that $v=...$ is just the expression for the value in brackets, it could be for scalars $v=\Psi(\vec{r}_i) \, \Psi(\vec{r}_j)$, vectors $v=\vec{a}(t_1) \cdot \vec{a}(t_2)$, etc Dec 3 '21 at 18:38