custom lattice commit in new branch

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)

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

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:

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