From a5518ef1ccd26242032432a6a874391c3dbcfe8f Mon Sep 17 00:00:00 2001 From: aaron4444 Date: Tue, 11 Jul 2017 11:09:52 -0400 Subject: [PATCH] custom lattice commit in new branch --- ScatterSim/CompositeNanoObjects.py | 173 +++++++++++++++++- ScatterSim/CustomLattices.py | 56 ++++++ ScatterSim/LatticeObjects.py | 6 +- ScatterSim/NanoObjects.py | 16 +- ScatterSim/old/MultiComponentModel.py | 248 +++++++++++++------------- 5 files changed, 363 insertions(+), 136 deletions(-) create mode 100644 ScatterSim/CustomLattices.py diff --git a/ScatterSim/CompositeNanoObjects.py b/ScatterSim/CompositeNanoObjects.py index 1fabf63..8929258 100644 --- a/ScatterSim/CompositeNanoObjects.py +++ b/ScatterSim/CompositeNanoObjects.py @@ -2,7 +2,7 @@ import numpy as np from copy import deepcopy # This file is where more complex nano objects can be stored -# These interface the same way as NanoObjects but since they're a little more complex +# These interface the same way as NanoObjects but since they're a little more complex # it makes sense to separate them from NanoObjects # TODO : Better simpler API for making composite object # like make_composite( objectslist, pargslist) @@ -233,7 +233,7 @@ def __init__(self, baseObject=None, linkerObject=None, pargs={}, seed=None): # you flip x or y from original shifts to move along edge axis # not a good explanation but some sort of personal bookkeeping for now... shiftfacs = [ - # top + # top [0,-1,1], [-1,0,1], [0,1,1], @@ -565,3 +565,172 @@ def __init__(self, pargs={}): objslist = [CylinderNanoObject, SphereNanoObject] super(CylinderBulgeNanoObject, self).__init__(objslist, pargslist, pargs=pargs) + + +class OctahedronColoredNanoObject(CompositeNanoObject): + """An octahedron object made of cylinders or like object. The canonical + (unrotated) version has the square cross-section in the x-y plane, with + corners pointing along +z and -z. The corners are on the x-axis and + y-axis. The edges are 45 degrees to the x and y axes. + The canonical (unrotated) version of the cylinders should be aligned along + the z axis. Replace cylinders with spheres, cylindrical shells, prolate + ellipsoids etc as you wish. + + It is best to just brute force define all the terms in one shot, + which I chose to do here. + notes about rotations: + - i need to dot the rotation matrix of the octahedron with each + individual cylinder's rotationo element + edgelength : length of edge + edgespread : how much to expand the rods by (not a good name) + positive is expansion, negative is compression + edgesep : separation of element from edge + rest of pargs are used for the cylinder object + ex : radius, height + linkerlength : if specified, adds linkers of specified length, + centered in between the octahedra + linkerradius : if linkerlength specified, will add linkers of this radius + rho_linker : if linkerlength specified, adds linkers of this density + (defaults to same density as cylinders in octahedra) + linkerobject : the object to use for the linkers (defaults to baseObject) + """ + def __init__(self, baseObject=None, linkerObject=None, pargs={}, seed=None): + if baseObject is None: + baseObject = CylinderNanoObject + if linkerObject is None: + linkerObject = baseObject + + # Set defaults + if 'edgeshift' not in pargs: + pargs['edgeshift'] = 0.0 + if 'edgespread' not in pargs: + pargs['edgespread'] = 0.0 + if 'linkerlength' in pargs: + addlinkers = True + else: + addlinkers = False + + # raise errors for undefined parameters + if 'edgelength' not in pargs: + raise ValueError("Need to specify an edgelength for this object") + + # these are slight shifts per cyl along the axis + # positive is away from COM and negative towards + shiftlabels = [ + # these correspond to the poslist + 'CYZ1', 'CXZ1', 'CYZ2', 'CXZ2', + 'CXY1', 'CXY4', 'CXY3', 'CXY2', + 'CYZ3', 'CXZ3', 'CYZ4', 'CXZ4', + 'linker1', 'linker2', 'linker3', 'linker4', + 'linker5', 'linker6', + ] + + # you flip x or y from original shifts to move along edge axis + # not a good explanation but some sort of personal bookkeeping for now... + shiftfacs = [ + # top + [0,-1,1], + [-1,0,1], + [0,1,1], + [1, 0,1], + # middle + [-1,1,0], + [-1,-1,0], + [1,-1,0], + [1,1,0], + # bottom + [0,1,-1], + [1, 0, -1], + [0,-1, -1], + [-1,0,-1] + ] + + for lbl1 in shiftlabels: + if lbl1 not in pargs: + pargs[lbl1] = 0. + + # calculate shift of COM from edgelength and edgespread + fac1 = np.sqrt(2)/2.*((.5*pargs['edgelength']) + pargs['edgespread']) + eL = pargs['edgelength'] + if addlinkers: + sL = pargs['linkerlength'] + + + poslist = [ + # eta, theta, phi, x0, y0, z0 + # top part + [0, 45, -90, 0, fac1, fac1], + [0, 45, 0, fac1, 0, fac1], + [0, 45, 90, 0, -fac1, fac1], + [0, -45, 0, -fac1, 0, fac1], + # now the flat part + [0, 90, 45, fac1, fac1, 0], + [0, 90, -45, fac1, -fac1, 0], + [0, 90, 45, -fac1, -fac1, 0], + [0, 90, -45, -fac1, fac1, 0], + # finally bottom part + [0, 45, -90, 0, -fac1,-fac1], + [0, 45, 0, -fac1, 0, -fac1], + [0, 45, 90, 0, fac1, -fac1], + [0, -45, 0, fac1, 0, -fac1], + ] + + if addlinkers: + poslist_linker = [ + # linkers + [0, 0, 0, 0, 0, eL/np.sqrt(2) + sL/2.], + [0, 0, 0, 0, 0, -eL/np.sqrt(2) - sL/2.], + [0, 90, 0, eL/np.sqrt(2) + sL/2.,0,0], + [0, 90, 0, -eL/np.sqrt(2) - sL/2.,0,0], + [0, 90, 90, 0, eL/np.sqrt(2) + sL/2., 0], + [0, 90, 90, 0, -eL/np.sqrt(2) - sL/2., 0], + ] + for row in poslist_linker: + poslist.append(row) + shiftfacs_linker = [ + [0,0,1], + [0,0,-1], + [1,0,0], + [-1,0,0], + [0,1,0], + [0,-1,0], + ] + for row in shiftfacs_linker: + shiftfacs.append(row) + + poslist = np.array(poslist) + shiftfacs = np.array(shiftfacs) + + # now add the shift factors + for i in range(len(poslist)): + poslist[i, 3:] += np.sqrt(2)/2.*shiftfacs[i]*pargs[shiftlabels[i]] + + + # need to create objslist and pargslist + objlist = list() + pargslist = list() + for i, pos in enumerate(poslist): + objlist.append(baseObject) + + eta, phi, theta, x0, y0, z0 = pos + pargstmp = dict() + pargstmp.update(pargs) + pargstmp['eta'] = eta + pargstmp['theta'] = theta + pargstmp['phi'] = phi + pargstmp['x0'] = x0 + pargstmp['y0'] = y0 + pargstmp['z0'] = z0 + + labeltmp = shiftlabels[i] + if 'linker' in labeltmp: + if 'rho_linker' in pargs: + pargstmp['rho_object'] = pargs['rho_linker'] + if 'linkerlength' in pargs: + pargstmp['height'] = pargs['linkerlength'] + if 'linkerradius' in pargs: + pargstmp['radius'] = pargs['linkerradius'] + + pargslist.append(pargstmp) + + super(OctahedronCylindersNanoObject, self).__init__(objlist, pargslist, pargs=pargs) diff --git a/ScatterSim/CustomLattices.py b/ScatterSim/CustomLattices.py new file mode 100644 index 0000000..2382e85 --- /dev/null +++ b/ScatterSim/CustomLattices.py @@ -0,0 +1,56 @@ +from .LatticeObjects import Lattice + + class BCCLatticeColor(Lattice): + def __init__(self, objects, lattice_spacing_a=1.0, sigma_D=0.01): + ''' cannot specify lattice_types or lattice_coordinates''' + # Define the lattice + symmetry = { + 'crystal family' : 'cubic', + 'crystal system' : 'cubic', + 'Bravais lattice' : 'I', + 'crystal class' : 'hexoctahedral', + 'point group' : 'm3m', + 'space group' : 'Im3m', + } + + lattice_coordinates = np.array([ (0.0, 0.0, 0.0), \ + (0.5, 0.5, 0.5), \ + ]) + lattice_types = [1,2] + positions = ['all'] + lattice_positions = ['corner', 'center'] + + + super(BCCLatticeColor, self).__init__(objects, lattice_spacing_a=lattice_spacing_a, + sigma_D=sigma_D, symmetry=symmetry, + lattice_positions=lattice_positions, + lattice_coordinates=lattice_coordinates, + lattice_types=lattice_types) + + def symmetry_factor(self, h, k, l): + """Returns the symmetry factor (0 for forbidden).""" + return 1 + if (h+k+l)%2==0: + return 2 + else: + return 0 + + def unit_cell_volume(self): + return self.lattice_spacing_a**3 + + def q_hkl(self, h, k, l): + """Determines the position in reciprocal space for the given reflection.""" + + prefactor = (2*np.pi/self.lattice_spacing_a) + qhkl_vector = np.array([ prefactor*h, \ + prefactor*k, \ + prefactor*l ]) + qhkl = np.sqrt( qhkl_vector[0]**2 + qhkl_vector[1]**2 + qhkl_vector[2]**2 ) + + return (qhkl, qhkl_vector) + + def q_hkl_length(self, h, k, l): + prefactor = (2*np.pi/self.lattice_spacing_a) + qhkl = prefactor*np.sqrt( h**2 + k**2 + l**2 ) + + return qhkl diff --git a/ScatterSim/LatticeObjects.py b/ScatterSim/LatticeObjects.py index 2341d92..92c8dc8 100644 --- a/ScatterSim/LatticeObjects.py +++ b/ScatterSim/LatticeObjects.py @@ -267,7 +267,7 @@ def sum_over_hkl(self, q, peak, max_hkl=6): fs = self.form_factor(qhkl_vector) term1 = (fs*fs.conjugate()).real - # if + # if #term2 = np.exp( -(self.sigma_D**2) * (qhkl**2) * (self.lattice_spacing_a**2) ) term2 = self.G_q(qhkl_vector[0], qhkl_vector[1], qhkl_vector[2]) term3 = peak( q-qhkl ) @@ -398,7 +398,7 @@ def to_string(self): s += " (a, b, c) = (%.3f,%.3f,%.3f) in nm\n" % (self.lattice_spacing_a,self.lattice_spacing_b,self.lattice_spacing_c) s += " (alpha, beta, gamma) = (%.3f,%.3f,%.3f) in radians\n" % (self.alpha,self.beta,self.gamma) - s += " = (%.2f,%.2f,%.2f) in degrees\n" % (degrees(self.alpha),degrees(self.beta),degrees(self.gamma)) + #s += " = (%.2f,%.2f,%.2f) in degrees\n" % (degrees(self.alpha),degrees(self.beta),degrees(self.gamma)) s += " volume = %.4f nm^3\n\n" % self.unit_cell_volume() s += " Objects:\n" for i, obj in enumerate(self.objects): @@ -432,7 +432,7 @@ def __init__(self, objects, lattice_spacing_a=1.0, sigma_D=0.01, lattice_coordin } if not isinstance(objects, list): objects = [objects] - + lattice_positions = ['placement']*len(objects) if lattice_coordinates is None: diff --git a/ScatterSim/NanoObjects.py b/ScatterSim/NanoObjects.py index 0a95609..8599570 100644 --- a/ScatterSim/NanoObjects.py +++ b/ScatterSim/NanoObjects.py @@ -866,7 +866,7 @@ class RandomizedNanoObject(NanoObject): whose radius follows the lognormal distribution, you would put: { 'radius' : { - 'distribution_type' : 'lognormal', + 'distribution_type' : 'lognormal', 'mean' : 1, # average val 'sigma' : .1, #std deviation }, @@ -895,7 +895,7 @@ def __init__(self, baseNanoObjectClass, pargs={}, argdict=None): NanoObject.__init__(self, pargs=pargs) if argdict is None: - argdict = dict(radius={'distribution' : 'gaussian', + argdict = dict(radius={'distribution' : 'gaussian', 'sigma' : .1, 'mean' : 1}) if 'distribution_num_points' not in pargs: pargs['distribution_num_points'] = 21 @@ -947,7 +947,7 @@ def build_objects(self): when summing run same object, rebuild and re-calculate since it's random, you don't want to cache either 'distribution_num_points' : 21, - {'radius' : {'distribution_type' : 'gaussian', + {'radius' : {'distribution_type' : 'gaussian', # parameters specific to distribution 'avg' : 1, # average val 'std' : .1, #std deviation @@ -964,13 +964,13 @@ def sample_point(self, dist_dict): ''' Sample a point from a distribution specified by dist_dict Currently supported: - 'distribution_type' (case insensitive) + 'distribution_type' (case insensitive) 'Gaussian' or 'normal' - parameters : + parameters : 'mean' : average of Gaussian (normal) distribution 'sigma' : standard deviation of Gaussian (normal) distribution - 'uniform' : - parameters : + 'uniform' : + parameters : 'low' : lower bound of distribution 'high' : upper bound of distribution 'lognormal' @@ -1253,7 +1253,7 @@ def form_factor(self, qvec): qr = np.hypot(qx, qy) - # NOTE : Numpy's sinc function adds + # NOTE : Numpy's sinc function adds # a factor of pi in we need to remove. # Why numpy... why??? >< F = 2*np.sinc(qz*H/2./np.pi)*j1(qr*R)/qr/R + 1j*0 diff --git a/ScatterSim/old/MultiComponentModel.py b/ScatterSim/old/MultiComponentModel.py index 6be8878..c5a3720 100644 --- a/ScatterSim/old/MultiComponentModel.py +++ b/ScatterSim/old/MultiComponentModel.py @@ -20,16 +20,17 @@ ################################################################### -from ScatterSim.BaseClasses import Potential, Model, arrays_equal -from ScatterSim.BaseClasses import radians, cos, sin +from ScatterSim.BaseClasses import Potential, Model, arrays_equal, value_at +from ScatterSim.BaseClasses import radians, cos, sin, degrees from ScatterSim import gamma import os, sys from cmath import exp as cexp +from numpy import sinc import random import numpy as np - +import matplotlib.pyplot as plt # Bessel functions of the first kind, orders 0 and 1 from scipy.special import j0, j1 @@ -492,31 +493,31 @@ def plot_form_factor_amplitude(self, qtuple, filename='form_factor_amplitude.png #q_zeros = np.zeros(len(q_list)) #int_list = self.form_factor_array(q_zeros,q_zeros,q_list) - pylab.rcParams['axes.labelsize'] = 30 - pylab.rcParams['xtick.labelsize'] = 'xx-large' - pylab.rcParams['ytick.labelsize'] = 'xx-large' - fig = pylab.figure() + plt.rcParams['axes.labelsize'] = 30 + plt.rcParams['xtick.labelsize'] = 'xx-large' + plt.rcParams['ytick.labelsize'] = 'xx-large' + fig = plt.figure() fig.subplots_adjust(left=0.14, bottom=0.15, right=0.94, top=0.94) - pylab.plot( q_list, int_list, color=(0,0,0), linewidth=3.0 ) + plt.plot( q_list, int_list, color=(0,0,0), linewidth=3.0 ) if ylog: - pylab.semilogy() + plt.semilogy() else: # Make y-axis scientific notation fig.gca().yaxis.major.formatter.set_scientific(True) fig.gca().yaxis.major.formatter.set_powerlimits((3,3)) - pylab.xlabel( r'$q \, (\mathrm{nm}^{-1})$' ) - pylab.ylabel( r'$F(q)$' ) + plt.xlabel( r'$q \, (\mathrm{nm}^{-1})$' ) + plt.ylabel( r'$F(q)$' ) - #xi, xf, yi, yf = pylab.axis() + #xi, xf, yi, yf = plt.axis() #yf = 5e5 #yi = 0 - #pylab.axis( [xi, xf, yi, yf] ) + #plt.axis( [xi, xf, yi, yf] ) - pylab.savefig( filename ) + plt.savefig( filename ) return int_list @@ -535,30 +536,30 @@ def plot_form_factor_intensity(self, qtuple, filename='form_factor_intensity.png int_list.append( self.form_factor_intensity(0,0,q) ) - pylab.rcParams['axes.labelsize'] = 30 - pylab.rcParams['xtick.labelsize'] = 'xx-large' - pylab.rcParams['ytick.labelsize'] = 'xx-large' - fig = pylab.figure() + plt.rcParams['axes.labelsize'] = 30 + plt.rcParams['xtick.labelsize'] = 'xx-large' + plt.rcParams['ytick.labelsize'] = 'xx-large' + fig = plt.figure() fig.subplots_adjust(left=0.14, bottom=0.15, right=0.94, top=0.94) - pylab.plot( q_list, int_list, color=(0,0,0), linewidth=3.0 ) + plt.plot( q_list, int_list, color=(0,0,0), linewidth=3.0 ) if ylog: - pylab.semilogy() + plt.semilogy() else: # Make y-axis scientific notation fig.gca().yaxis.major.formatter.set_scientific(True) fig.gca().yaxis.major.formatter.set_powerlimits((3,3)) - pylab.xlabel( r'$q \, (\mathrm{nm}^{-1})$' ) - pylab.ylabel( r'$F(q)$' ) + plt.xlabel( r'$q \, (\mathrm{nm}^{-1})$' ) + plt.ylabel( r'$F(q)$' ) - #xi, xf, yi, yf = pylab.axis() + #xi, xf, yi, yf = plt.axis() #yf = 5e5 - #pylab.axis( [xi, xf, yi, yf] ) + #plt.axis( [xi, xf, yi, yf] ) - pylab.savefig( filename ) + plt.savefig( filename ) return int_list @@ -1837,7 +1838,7 @@ class PolydisperseRodNanoObject(NanoObject): experimental. To average over orientations, a new rotation matrix will need to be - computed. + computed. """ @@ -3178,7 +3179,8 @@ def V(self, in_x, in_y, in_z, rotation_elements=None): def form_factor(self, qx, qy, qz): """Returns the complex-amplitude of the form factor at the given q-coordinates.""" - + import numpy + from numpy import sinc qx, qy, qz = self.rotate_coord(qx, qy, qz) R = self.pargs['radius'] @@ -5476,7 +5478,7 @@ def __init__(self, baseObject=None, pargs={}, seed=None): # you flip x or y from original shifts to move along edge axis # not a good explanation but some sort of personal bookkeeping for now... shiftfacs = np.array([ - # top + # top [0,-1,1], [-1,0,1], [0,1,1], @@ -5708,7 +5710,7 @@ def __init__(self, baseObject=None, pargs={}, seed=None): # you flip x or y from original shifts to move along edge axis # not a good explanation but some sort of personal bookkeeping for now... shiftfacs = np.array([ - # top + # top [0,-1,1], [-1,0,1], [0,1,1], @@ -6637,21 +6639,21 @@ def plot(self, plot_width=1.0, num_points=200, filename='peak.png', ylog=False): q_list = np.linspace( -plot_width, plot_width, num_points ) int_list = self.val_array( q_list, 0.0 ) - pylab.rcParams['axes.labelsize'] = 30 - pylab.rcParams['xtick.labelsize'] = 'xx-large' - pylab.rcParams['ytick.labelsize'] = 'xx-large' - fig = pylab.figure() + plt.rcParams['axes.labelsize'] = 30 + plt.rcParams['xtick.labelsize'] = 'xx-large' + plt.rcParams['ytick.labelsize'] = 'xx-large' + fig = plt.figure() fig.subplots_adjust(left=0.14, bottom=0.15, right=0.94, top=0.94) - pylab.plot( q_list, int_list, color=(0,0,0), linewidth=3.0 ) + plt.plot( q_list, int_list, color=(0,0,0), linewidth=3.0 ) if ylog: - pylab.semilogy() + plt.semilogy() - pylab.xlabel( r'$q \, (\mathrm{nm}^{-1})$', size=30 ) - pylab.ylabel( 'Intensity (a.u.)', size=30 ) + plt.xlabel( r'$q \, (\mathrm{nm}^{-1})$', size=30 ) + plt.ylabel( 'Intensity (a.u.)', size=30 ) - pylab.savefig( filename ) + plt.savefig( filename ) return int_list @@ -6909,7 +6911,7 @@ def iterate_over_hkl_1d(self, max_hkl=6, tolerance=1e-10): # Search dictionary, see if this qhkl already exists found = False - for q in sorted(hkl_list_1d.iterkeys()): + for q in sorted(hkl_list_1d.keys()): if abs(qhkl-q)