Shifting ExampleΒΆ

This notebook shows how to use the shifting option on PMFTXYZ to get high resolution views of PMFT features that are not centered.

[1]:
import numpy as np
from freud import box, pmft

from scipy.interpolate import griddata
from scipy.interpolate import RegularGridInterpolator

import warnings
warnings.simplefilter('ignore')

import matplotlib.pyplot as plt

First we load in our data. The particles used here are implemented with a simple Weeks-Chandler-Andersen isotropic pair potential, so particle orientation is not meaningful.

[2]:
pos_data = np.load('data/XYZ/positions.npy').astype(np.float32)
box_data = np.load('data/XYZ/boxes.npy').astype(np.float32)

We calculate the PMFT the same way as shown in other examples first

[3]:
window = 2**(1/6) # The size of the pmft calculation

res = (100,100,100)
pmft_arr = np.zeros(res)

mypmft = pmft.PMFTXYZ(x_max=window, y_max=window, z_max=window,
                      n_x=res[0], n_y=res[1], n_z=res[2])

# This data is for isotropic particles, so we will just make some unit quaternions
# to use as the orientations
quats = np.zeros((pos_data.shape[1],4)).astype(np.float32)
quats[:,0] = 1

for i in range(10, pos_data.shape[0]):
    l_box = box_data[i]
    l_pos = pos_data[i]
    mypmft.accumulate(l_box, l_pos, quats, l_pos, quats)

unshifted = np.copy(mypmft.PMFT)

x = mypmft.X
y = mypmft.Y
z = mypmft.Z

When we plot a centered slice of the XYZ pmft, we see that a number of wells are present at some distance from the origin

[4]:
%matplotlib inline

plt.figure(figsize=(10,10))
plt.imshow(unshifted[int(res[0]/2),:,:])
plt.colorbar()
plt.show()
../../_images/examples_module_intros_PMFT-PMFTXYZ_Shift_Example_7_0.png

If we want a closer look at the details of those wells, then we could increase the PMFT resolution. But this will increase the computational cost by a lot, and we are wasting a big percentage of the pixels.

This use case is why the shiftvec argument was implemented. Now we will do the same calculation, but we will use a much smaller window centered on on of the wells.

To do this we need to pass a vector into the PMFTXYZ construction. The window will be centered on this vector.

[5]:
shiftvec = [0.82, 0.82, 0]

window = 2**(1/6)/6 # Smaller window for the shifted case

res = (50,50,50)
pmft_arr = np.zeros(res)

mypmft = pmft.PMFTXYZ(x_max=window, y_max=window, z_max=window,
                      n_x=res[0], n_y=res[1], n_z=res[2], shiftvec=shiftvec)

# This data is for isotropic particles, so we will just make some unit quaternions
# to use as the orientations
quats = np.zeros((pos_data.shape[1],4)).astype(np.float32)
quats[:,0] = 1

for i in range(10,pos_data.shape[0]):
    l_box = box_data[i]
    l_pos = pos_data[i]
    mypmft.accumulate(l_box, l_pos, quats, l_pos, quats)

shifted = np.copy(mypmft.PMFT)

x = mypmft.X
y = mypmft.Y
z = mypmft.Z

Now the PMFT is a high resolution close up of one of the bonding wells. Note that as you increase the sampling resolution, you need to increase your number of samples because there is less averaging in each bin

[6]:
%matplotlib inline

plt.figure(figsize=(10,10))
plt.imshow(shifted[int(res[2]/2),:,:])
plt.colorbar()
plt.show()
../../_images/examples_module_intros_PMFT-PMFTXYZ_Shift_Example_11_0.png