Changeset 2a39ca4 in sasmodels

Timestamp:

Sep 25, 2018 12:20:12 PM (6 years ago)

Author:

Paul Kienzle <pkienzle@…>

Branches:

ticket_1156

Children:

5778279

Parents:

642046e

Message:

sample from sc/bcc/fcc lattice in real space calculation

File:

• explore/realspace.py (modified) (13 diffs)

explore/realspace.py

@@ -78 +78 @@
 
 
+I3 = np.matrix([[1., 0, 0], [0, 1, 0], [0, 0, 1]])
+
 class Shape:
-    rotation = np.matrix([[1., 0, 0], [0, 1, 0], [0, 0, 1]])
+    rotation = I3
     center = np.array([0., 0., 0.])[:, None]
     r_max = None
@@ -96 +98 @@
 
     def rotate(self, theta, phi, psi):
-        self.rotation = rotation(theta, phi, psi) * self.rotation
+        if theta != 0. or phi != 0. or psi != 0.:
+            self.rotation = rotation(theta, phi, psi) * self.rotation
         return self
 
@@ -104 +107 @@
 
     def _adjust(self, points):
-        points = np.asarray(self.rotation * np.matrix(points.T)) + self.center
+        if self.rotation is I3:
+            points = points.T + self.center
+        else:
+            points = np.asarray(self.rotation * np.matrix(points.T)) + self.center
         return points.T
 
@@ -669 +675 @@
     Iq = 100 * np.ones_like(qx)
     data = Data2D(x=qx, y=qy, z=Iq, dx=None, dy=None, dz=np.sqrt(Iq))
-    data.x_bins = qx[0,:]
-    data.y_bins = qy[:,0]
+    data.x_bins = qx[0, :]
+    data.y_bins = qy[:, 0]
     data.filename = "fake data"
 
@@ -695 +701 @@
     return shape, fn, fn_xy
 
-def build_sphere(radius=125, rho=2):
+DEFAULT_SPHERE_RADIUS = 125
+DEFAULT_SPHERE_CONTRAST = 2
+def build_sphere(radius=DEFAULT_SPHERE_RADIUS, rho=DEFAULT_SPHERE_CONTRAST):
     shape = TriaxialEllipsoid(radius, radius, radius, rho)
     fn = lambda q: sphere_Iq(q, radius)*rho**2
@@ -751 +759 @@
     return shape, fn, fn_xy
 
-def build_cubic_lattice(shape, nx=1, ny=1, nz=1, dx=2, dy=2, dz=2,
-                  shuffle=0, rotate=0):
+def build_sc_lattice(shape, nx=1, ny=1, nz=1, dx=2, dy=2, dz=2,
+                        shuffle=0, rotate=0):
     a, b, c = shape.dims
-    shapes = [copy(shape)
+    corners= [copy(shape)
               .shift((ix+(randn() if shuffle < 0.3 else rand())*shuffle)*dx*a,
                      (iy+(randn() if shuffle < 0.3 else rand())*shuffle)*dy*b,
@@ -762 +770 @@
               for iy in range(ny)
               for iz in range(nz)]
-    lattice = Composite(shapes)
+    lattice = Composite(corners)
     return lattice
 
+def build_bcc_lattice(shape, nx=1, ny=1, nz=1, dx=2, dy=2, dz=2,
+                      shuffle=0, rotate=0):
+    a, b, c = shape.dims
+    corners = [copy(shape)
+               .shift((ix+(randn() if shuffle < 0.3 else rand())*shuffle)*dx*a,
+                      (iy+(randn() if shuffle < 0.3 else rand())*shuffle)*dy*b,
+                      (iz+(randn() if shuffle < 0.3 else rand())*shuffle)*dz*c)
+               .rotate(*((randn(3) if rotate < 30 else rand(3))*rotate))
+               for ix in range(nx)
+               for iy in range(ny)
+               for iz in range(nz)]
+    centers = [copy(shape)
+               .shift((ix+0.5+(randn() if shuffle < 0.3 else rand())*shuffle)*dx*a,
+                      (iy+0.5+(randn() if shuffle < 0.3 else rand())*shuffle)*dy*b,
+                      (iz+0.5+(randn() if shuffle < 0.3 else rand())*shuffle)*dz*c)
+               .rotate(*((randn(3) if rotate < 30 else rand(3))*rotate))
+               for ix in range(nx)
+               for iy in range(ny)
+               for iz in range(nz)]
+    lattice = Composite(corners + centers)
+    return lattice
+
+def build_fcc_lattice(shape, nx=1, ny=1, nz=1, dx=2, dy=2, dz=2,
+                      shuffle=0, rotate=0):
+    a, b, c = shape.dims
+    corners = [copy(shape)
+               .shift((ix+(randn() if shuffle < 0.3 else rand())*shuffle)*dx*a,
+                      (iy+(randn() if shuffle < 0.3 else rand())*shuffle)*dy*b,
+                      (iz+(randn() if shuffle < 0.3 else rand())*shuffle)*dz*c)
+               .rotate(*((randn(3) if rotate < 30 else rand(3))*rotate))
+               for ix in range(nx)
+               for iy in range(ny)
+               for iz in range(nz)]
+    faces_a = [copy(shape)
+               .shift((ix+0.0+(randn() if shuffle < 0.3 else rand())*shuffle)*dx*a,
+                      (iy+0.5+(randn() if shuffle < 0.3 else rand())*shuffle)*dy*b,
+                      (iz+0.5+(randn() if shuffle < 0.3 else rand())*shuffle)*dz*c)
+               .rotate(*((randn(3) if rotate < 30 else rand(3))*rotate))
+               for ix in range(nx)
+               for iy in range(ny)
+               for iz in range(nz)]
+    faces_b = [copy(shape)
+               .shift((ix+0.5+(randn() if shuffle < 0.3 else rand())*shuffle)*dx*a,
+                      (iy+0.0+(randn() if shuffle < 0.3 else rand())*shuffle)*dy*b,
+                      (iz+0.5+(randn() if shuffle < 0.3 else rand())*shuffle)*dz*c)
+               .rotate(*((randn(3) if rotate < 30 else rand(3))*rotate))
+               for ix in range(nx)
+               for iy in range(ny)
+               for iz in range(nz)]
+    faces_c = [copy(shape)
+               .shift((ix+0.5+(randn() if shuffle < 0.3 else rand())*shuffle)*dx*a,
+                      (iy+0.5+(randn() if shuffle < 0.3 else rand())*shuffle)*dy*b,
+                      (iz+0.0+(randn() if shuffle < 0.3 else rand())*shuffle)*dz*c)
+               .rotate(*((randn(3) if rotate < 30 else rand(3))*rotate))
+               for ix in range(nx)
+               for iy in range(ny)
+               for iz in range(nz)]
+    lattice = Composite(corners + faces_a + faces_b + faces_c)
+    return lattice
 
 SHAPE_FUNCTIONS = OrderedDict([
@@ -775 +842 @@
 ])
 SHAPES = list(SHAPE_FUNCTIONS.keys())
+LATTICE_FUNCTIONS = OrderedDict([
+    ("sc", build_sc_lattice),
+    ("bcc", build_bcc_lattice),
+    ("fcc", build_fcc_lattice),
+])
+LATTICE_TYPES = list(LATTICE_FUNCTIONS.keys())
 
 def check_shape(title, shape, fn=None, show_points=False,
@@ -783 +856 @@
     r = shape.r_bins(q, r_step=r_step)
     sampling_density = samples / shape.volume
+    print("sampling points")
     rho, points = shape.sample(sampling_density)
+    print("calculating Pr")
     t0 = time.time()
     Pr = calc_Pr(r, rho-rho_solvent, points)
@@ -806 +881 @@
     Qx, Qy = np.meshgrid(qx, qy)
     sampling_density = samples / shape.volume
+    print("sampling points")
     t0 = time.time()
     rho, points = shape.sample(sampling_density)
@@ -844 +920 @@
                         help='lattice size')
     parser.add_argument('-z', '--spacing', type=str, default='2,2,2',
-                        help='lattice spacing')
+                        help='lattice spacing (relative to shape)')
+    parser.add_argument('-t', '--type', choices=LATTICE_TYPES,
+                        default=LATTICE_TYPES[0],
+                        help='lattice type')
     parser.add_argument('-r', '--rotate', type=float, default=0.,
                         help="rotation relative to lattice, gaussian < 30 degrees, uniform otherwise")
@@ -860 +939 @@
     shuffle, rotate = opts.shuffle, opts.rotate
     shape, fn, fn_xy = SHAPE_FUNCTIONS[opts.shape](**pars)
+    view = tuple(float(v) for v in opts.view.split(','))
     if nx > 1 or ny > 1 or nz > 1:
-        shape = build_cubic_lattice(shape, nx, ny, nz, dx, dy, dz, shuffle, rotate)
+        print("building %s lattice"%opts.type)
+        lattice = LATTICE_FUNCTIONS[opts.type]
+        shape = lattice(shape, nx, ny, nz, dx, dy, dz, shuffle, rotate)
+        # If comparing a sphere in a cubic lattice, compare against the
+        # corresponding paracrystalline model.
+        if opts.shape == "sphere" and dx == dy == dz:
+            radius = pars.get('radius', DEFAULT_SPHERE_RADIUS)
+            model_name = opts.type + "_paracrystal"
+            model_pars = {
+                "scale": 1.,
+                "background": 0.,
+                "lattice_spacing": 2*radius*dx,
+                "lattice_distortion": shuffle,
+                "radius": radius,
+                "sld": pars.get('rho', DEFAULT_SPHERE_CONTRAST),
+                "sld_solvent": 0.,
+                "theta": view[0],
+                "phi": view[1],
+                "psi": view[2],
+            }
+            fn, fn_xy = wrap_sasmodel(model_name, **model_pars)
     title = "%s(%s)" % (opts.shape, " ".join(opts.pars))
     if opts.dim == 1:
@@ -867 +967 @@
                     mesh=opts.mesh, qmax=opts.qmax, samples=opts.samples)
     else:
-        view = tuple(float(v) for v in opts.view.split(','))
         check_shape_2d(title, shape, fn_xy, view=view, show_points=opts.plot,
                        mesh=opts.mesh, qmax=opts.qmax, samples=opts.samples)