<jupyter_start><jupyter_text>This code uses the freud package, alongwith PyLammps, to calculate the averaged structure factor over a user-defined range of snapshots in a lammps trajectory file. TODO: Implement structure factor instead of RDF<jupyter_code>import numpy as np
import time
from cycler import cycler
import matplotlib
from matplotlib import pyplot
from matplotlib.colors import colorConverter
myHexColors = ["#30A2DA", "#FC4F30", "#E5AE38", "#6D904F", "#9757DB", "#188487", "#FF7F00", "#9A2C66", "#626DDA", "#8B8B8B"]
myColorCycle = [colorConverter.to_rgb(i) for i in myHexColors]
from freud import parallel
parallel.setNumThreads(4)

from ipywidgets import IntProgress
from IPython.display import display

# matplotlib.rcParams.update({'font.sans-serif': 'Helvetica'})
# matplotlib.rcParams.update({'font.family': 'sans-serif'})
matplotlib.rcParams.update({'font.size': 24})
matplotlib.rcParams.update({'axes.labelsize': 24})
matplotlib.rcParams.update({'xtick.labelsize': 18})
matplotlib.rcParams.update({'ytick.labelsize': 18})
matplotlib.rcParams.update({'savefig.pad_inches': 0.025})
matplotlib.rcParams.update({"axes.prop_cycle": cycler(color=myColorCycle)})
matplotlib.rcParams.update({"lines.linewidth": 2})

# lammps
from lammps import lammps, PyLammps
lmp = lammps()
L = PyLammps(ptr=lmp) # Get PyLammps<jupyter_output><empty_output><jupyter_text>All the user-defined stuff is in the next cell. Update the name of your trajectory file accordingly. It should be inside the subfolder traj. If you specify a frame that is not present in the lammps trajectory file, PyLammps will throw an error. TODO: Read these from an input file?<jupyter_code># User-defined stuff
system = 'liq'
typeID = 1
frameNo = 500000 #500000
endFrame = 515000 #1000000
dumpFreq = 5000 #1000 #5000
nFrames = int((endFrame - frameNo)/dumpFreq) + 1
# Dumpfile name to read from 
dumpfile = 'traj/dump-' + system + '.lammpstrj'
# Define the no. of grid points in each dim: 2*g+1
g = 10<jupyter_output><empty_output><jupyter_text>Average the RDF over the frames specified by frameNo and endFrame (inclusive of both). If you specify a frame that is not present in the lammps trajectory file, PyLammps will throw an error.If you run the cell below more than once, make sure to uncomment the PyLammps command for reading in the data file, or else you will get a Lammps error.<jupyter_code># import the freud object
# the rdf module is contained in the density module
from freud import box, density

# Current frame
frame = frameNo

# compute the rdf for for all frames except the first (your syntax will vary based on your reader)
myProgressBar = IntProgress(min=1,max=nFrames)
display(myProgressBar)

# Lammps init
# load the data
# Read in the lammps input file
lmp.file('lammpsFiles/inp.'+ system) # LAMMPS input file for the system
# Get the number of atoms of one typeID
natoms = 0
L.read_dump(dumpfile, frameNo, "x y z box yes replace yes")
for i in range(L.atoms.natoms):
    if L.atoms[i].type==typeID:
        natoms += 1
# Define numPy coordinate array; here natoms is the no. of atoms of type 1
x = np.zeros((natoms,3), dtype=float)

# Get the simulation box from PyLammps
box_limits = [L.system.xhi-L.system.xlo, L.system.yhi-L.system.ylo, L.system.zhi-L.system.zlo] # box from lammps
# create the freud box object
fbox = box.Box(Lx=box_limits[0], Ly=box_limits[1], Lz=box_limits[2]) 

# create the structure factor object
# g defines the no. of grid points = 2*g+1
strucFactor = density.RDF(box=fbox, g) 

start_time = time.time()
for i in range(nFrames):
    # Progress bar
    myProgressBar.value = int(frame/dumpFreq)
    # --------------------------
    # Lammps stuff
    # read box, position data from the lammps dumpfile
    L.read_dump(dumpfile, frame, "x y z box yes replace yes")
    # Get coord
    # Assigning values to numpy array x
    for i in range(L.atoms.natoms):
        if L.atoms[i].type==typeID:
            for k in range(3):
                x[i][k] = L.atoms[i].position[k]          
    box_limits = [L.system.xhi-L.system.xlo, L.system.yhi-L.system.ylo, L.system.zhi-L.system.zlo] # box from lammps
    # --------------------------
    # create the freud box object
    fbox = box.Box(Lx=box_limits[0], Ly=box_limits[1], Lz=box_limits[2])    

    # accumulate
    rdf.accumulate(fbox, x, x)

    # ----------------------------
    # Update frame number
    print("Frame number ", frame)
    frame = frame + dumpFreq
stop_time = time.time()
print("time to calc {} frames = {}".format(nFrames, stop_time-start_time))
print("speed of calc: {} (frames/sec)".format((nFrames)/(stop_time-start_time)))

# get the center of the histogram bins
r = rdf.getR()
# get the value of the histogram bins
y = rdf.getRDF()

pyplot.figure(figsize=(8, 5))
pyplot.title("RDF - Avg Frames")
pyplot.plot(r, y)
pyplot.xlabel(r"$r$")
pyplot.ylabel(r"$g\left(r\right)$")
pyplot.tight_layout()
pyplot.show()
pyplot.savefig('averaged_RDF.png')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12912308930008046, 'fraction_numerical': 0.021922767497988736, 'mean_word_length': 4.158713692946058, 'pattern_counts': {'":': 2, '<': 13, '<?xml version=': 0, '>': 13, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '251315', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 1797, 'n_tokens_neox': 1703, 'n_words': 526, 'path': '.ipynb_checkpoints/structureFactor-checkpoint.ipynb', 'repo_name': 'amritagos/structureFactor', 'size': '4972'}