Add Orientation Correlation Diagram

[Eric Harper]

Evaluate current bond order code, then add whatever the orientation correlation diagram is (may also want to look at the correlation function code to see if that helps at all)

[Eric Harper]

(Comments don’t line up with master right now, some were old and I added some to dos)

This starts at line 181 in BondOrder.cc:

// This needs to be moved to where the parallel_for is (look at pmft code to see example, will need to make sure that variable names are changed/correct)

// r is the parallel index
void operator()( const blocked_range<size_t>& r ) const
            {
            // we pre-calculate the inverses since division is more expensive
            float dt_inv = 1.0f / m_dt;
            float dp_inv = 1.0f / m_dp;
            float rmaxsq = m_rmax * m_rmax;
            // create object which indexes into the theta, phi histogram
            Index2D sa_i = Index2D(m_nbins_t, m_nbins_p);

            // this code creates per-thread memory, won't need to edit
            bool exists;
            m_local_bin_counts.local(exists);
            if (! exists)
                {
                m_local_bin_counts.local() = new unsigned int [m_nbins_t*m_nbins_p];
                memset((void*)m_local_bin_counts.local(), 0, sizeof(unsigned int)*m_nbins_t*m_nbins_p);
                }

        // This is the actual for loop    
        for(size_t i=r.begin(); i!=r.end(); ++i)
                {

                // get the position and orientation of particle i
                // do here because it's faster
                vec3<float> ref_pos = m_ref_points[i];
                quat<float> ref_q(m_ref_orientations[i]);

                //loop over neighbors
                locality::NearestNeighbors::iteratorneighbor it = m_nn->iterneighbor(i);
                for (unsigned int j = it.begin(); !it.atEnd(); j = it.next())
                    {

                    //compute r between the two particles
                    vec3<float> delta = m_box.wrap(m_points[j] - ref_pos);

                    float rsq = dot(delta, delta);
                // this is a "hack" that prevents particle i from binning itself    
                if (rsq > 1e-6)
                        {
                        // you will need to add this in 
                        // quat<float> q(m_orientations[j]);
                        // this should probably also be
                        // vec3<float> v = rotate(conj(ref_q), delta);
                        // you also need to add
                        // v = rotate(q, v);
                        vec3<float> v(delta);
                        v = rotate(conj(ref_q), v);
                        // get theta, phi
                        float theta = atan2f(v.y, v.x);
                        // these two lines need to be updated to the method in pmftR12
                        theta = (theta < 0) ? theta+2*M_PI : theta;
                        theta = (theta > 2*M_PI) ? theta-2*M_PI : theta;
                        float phi = atan2f(sqrt(v.x*v.x + v.y*v.y), v.z);
                        phi = (phi < 0) ? phi+2*M_PI : phi;
                        phi = (phi > 2*M_PI) ? phi-2*M_PI : phi;
                        // bin the point
                        float bint = floorf(theta * dt_inv);
                        float binp = floorf(phi * dp_inv);
                        // fast float to int conversion with truncation
                        #ifdef __SSE2__
                        unsigned int ibint = _mm_cvtt_ss2si(_mm_load_ss(&bint));
                        unsigned int ibinp = _mm_cvtt_ss2si(_mm_load_ss(&binp));
                        #else
                        unsigned int ibint = (unsigned int)(bint);
                        unsigned int ibinp = (unsigned int)(binp);
                        #endif

                        // increment the bin
                        if ((ibint < m_nbins_t) && (ibinp < m_nbins_p))
                            {
                            ++m_local_bin_counts.local()[sa_i(ibint, ibinp)];
                            }
                        }
                    }
                }
            }

[Eric Harper]

Here is what a script would look like

from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from matplotlib.colormap import cm
import freud

freud_handle = freud.order.BondOrder(rmax,
                                                              k,
                                                              n,
                                                              n_bins_theta,
                                                              n_bins_phi)
for frame in frames:
    freud_handle.accumulate(frame.box,
                                              frame.position,
                                              frame.orientation,
                                              frame.position,
                                              frame.orientation)
bond_order_array = self.compute.getBondOrder()
theta_array = self.compute.getTheta()
phi_array = self.compute.getPhi()
bond_order_array /= np.nanmax(bond_order_array)

# need to remap variables to plot in matplotlib
map_arr = np.zeros(shape=(bond_order_array.shape[1], bond_order_array.shape[0]), dtype=np.float32)
for i in range(bond_order_array.shape[1]):
    for j in range(bond_order_array.shape[0]):
    # this might actually be a transpose
    map_array[i,j] = bond_order_array[j,i]
theta = np.linspace(0.0, 2.0 * np.pi, n_bins_t+1)
phi = np.linspace(0.0, np.pi, n_bins_p+1)
x = np.outer(np.cos(theta), np.sin(phi))
y = np.outer(np.sin(theta), np.sin(phi))
z = np.outer(np.ones_like(theta), np.cos(phi))
# I'm not entirely sure if this code works
fig = plt.figure()
ax = fig.gca(projection='3d')
surf = ax.plot_surface(x, y, z, rstride=1, cstride=1, facecolors=cm.hot(map_array))
# ax.axis("off")
ax.grid(False)
for a in (ax.w_xaxis, ax.w_yaxis, ax.w_zaxis):
    for t in a.get_ticklines()+a.get_ticklabels():
            t.set_visible(False)
        a.line.set_visible(False)
        a.pane.set_visible(False)

[Eric Harper]

@erteich @djulia did this ever happen? Was part of hackathon, I don't recall if this was resolved.

[Bradley Dice]

@erteich @djulia Just following up to see the status of this, and whether there is work to be done or if the issue can be closed.

[Erin Teich]

Julia and I wrote these tools in the branch "ocd" during a Hackathon and that pull request was rejected by Eric in June 2016, who said he would merge it in manually:

https://bitbucket.org/glotzer/freud/pull-requests/67/ocd/diff

[Bradley Dice]

@erteich Thanks so much!

[Bradley Dice]

@vramasub confirmed that the changes were made manually as described.