Changeset 2a39ca4 in sasmodels


Ignore:
Timestamp:
Sep 25, 2018 12:20:12 PM (3 months ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
ticket_1156
Children:
5778279
Parents:
642046e
Message:

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

File:
1 edited

Legend:

Unmodified
Added
Removed
  • explore/realspace.py

    r362a66f r2a39ca4  
    7878 
    7979 
     80I3 = np.matrix([[1., 0, 0], [0, 1, 0], [0, 0, 1]]) 
     81 
    8082class Shape: 
    81     rotation = np.matrix([[1., 0, 0], [0, 1, 0], [0, 0, 1]]) 
     83    rotation = I3 
    8284    center = np.array([0., 0., 0.])[:, None] 
    8385    r_max = None 
     
    9698 
    9799    def rotate(self, theta, phi, psi): 
    98         self.rotation = rotation(theta, phi, psi) * self.rotation 
     100        if theta != 0. or phi != 0. or psi != 0.: 
     101            self.rotation = rotation(theta, phi, psi) * self.rotation 
    99102        return self 
    100103 
     
    104107 
    105108    def _adjust(self, points): 
    106         points = np.asarray(self.rotation * np.matrix(points.T)) + self.center 
     109        if self.rotation is I3: 
     110            points = points.T + self.center 
     111        else: 
     112            points = np.asarray(self.rotation * np.matrix(points.T)) + self.center 
    107113        return points.T 
    108114 
     
    669675    Iq = 100 * np.ones_like(qx) 
    670676    data = Data2D(x=qx, y=qy, z=Iq, dx=None, dy=None, dz=np.sqrt(Iq)) 
    671     data.x_bins = qx[0,:] 
    672     data.y_bins = qy[:,0] 
     677    data.x_bins = qx[0, :] 
     678    data.y_bins = qy[:, 0] 
    673679    data.filename = "fake data" 
    674680 
     
    695701    return shape, fn, fn_xy 
    696702 
    697 def build_sphere(radius=125, rho=2): 
     703DEFAULT_SPHERE_RADIUS = 125 
     704DEFAULT_SPHERE_CONTRAST = 2 
     705def build_sphere(radius=DEFAULT_SPHERE_RADIUS, rho=DEFAULT_SPHERE_CONTRAST): 
    698706    shape = TriaxialEllipsoid(radius, radius, radius, rho) 
    699707    fn = lambda q: sphere_Iq(q, radius)*rho**2 
     
    751759    return shape, fn, fn_xy 
    752760 
    753 def build_cubic_lattice(shape, nx=1, ny=1, nz=1, dx=2, dy=2, dz=2, 
    754                   shuffle=0, rotate=0): 
     761def build_sc_lattice(shape, nx=1, ny=1, nz=1, dx=2, dy=2, dz=2, 
     762                        shuffle=0, rotate=0): 
    755763    a, b, c = shape.dims 
    756     shapes = [copy(shape) 
     764    corners= [copy(shape) 
    757765              .shift((ix+(randn() if shuffle < 0.3 else rand())*shuffle)*dx*a, 
    758766                     (iy+(randn() if shuffle < 0.3 else rand())*shuffle)*dy*b, 
     
    762770              for iy in range(ny) 
    763771              for iz in range(nz)] 
    764     lattice = Composite(shapes) 
     772    lattice = Composite(corners) 
    765773    return lattice 
    766774 
     775def build_bcc_lattice(shape, nx=1, ny=1, nz=1, dx=2, dy=2, dz=2, 
     776                      shuffle=0, rotate=0): 
     777    a, b, c = shape.dims 
     778    corners = [copy(shape) 
     779               .shift((ix+(randn() if shuffle < 0.3 else rand())*shuffle)*dx*a, 
     780                      (iy+(randn() if shuffle < 0.3 else rand())*shuffle)*dy*b, 
     781                      (iz+(randn() if shuffle < 0.3 else rand())*shuffle)*dz*c) 
     782               .rotate(*((randn(3) if rotate < 30 else rand(3))*rotate)) 
     783               for ix in range(nx) 
     784               for iy in range(ny) 
     785               for iz in range(nz)] 
     786    centers = [copy(shape) 
     787               .shift((ix+0.5+(randn() if shuffle < 0.3 else rand())*shuffle)*dx*a, 
     788                      (iy+0.5+(randn() if shuffle < 0.3 else rand())*shuffle)*dy*b, 
     789                      (iz+0.5+(randn() if shuffle < 0.3 else rand())*shuffle)*dz*c) 
     790               .rotate(*((randn(3) if rotate < 30 else rand(3))*rotate)) 
     791               for ix in range(nx) 
     792               for iy in range(ny) 
     793               for iz in range(nz)] 
     794    lattice = Composite(corners + centers) 
     795    return lattice 
     796 
     797def build_fcc_lattice(shape, nx=1, ny=1, nz=1, dx=2, dy=2, dz=2, 
     798                      shuffle=0, rotate=0): 
     799    a, b, c = shape.dims 
     800    corners = [copy(shape) 
     801               .shift((ix+(randn() if shuffle < 0.3 else rand())*shuffle)*dx*a, 
     802                      (iy+(randn() if shuffle < 0.3 else rand())*shuffle)*dy*b, 
     803                      (iz+(randn() if shuffle < 0.3 else rand())*shuffle)*dz*c) 
     804               .rotate(*((randn(3) if rotate < 30 else rand(3))*rotate)) 
     805               for ix in range(nx) 
     806               for iy in range(ny) 
     807               for iz in range(nz)] 
     808    faces_a = [copy(shape) 
     809               .shift((ix+0.0+(randn() if shuffle < 0.3 else rand())*shuffle)*dx*a, 
     810                      (iy+0.5+(randn() if shuffle < 0.3 else rand())*shuffle)*dy*b, 
     811                      (iz+0.5+(randn() if shuffle < 0.3 else rand())*shuffle)*dz*c) 
     812               .rotate(*((randn(3) if rotate < 30 else rand(3))*rotate)) 
     813               for ix in range(nx) 
     814               for iy in range(ny) 
     815               for iz in range(nz)] 
     816    faces_b = [copy(shape) 
     817               .shift((ix+0.5+(randn() if shuffle < 0.3 else rand())*shuffle)*dx*a, 
     818                      (iy+0.0+(randn() if shuffle < 0.3 else rand())*shuffle)*dy*b, 
     819                      (iz+0.5+(randn() if shuffle < 0.3 else rand())*shuffle)*dz*c) 
     820               .rotate(*((randn(3) if rotate < 30 else rand(3))*rotate)) 
     821               for ix in range(nx) 
     822               for iy in range(ny) 
     823               for iz in range(nz)] 
     824    faces_c = [copy(shape) 
     825               .shift((ix+0.5+(randn() if shuffle < 0.3 else rand())*shuffle)*dx*a, 
     826                      (iy+0.5+(randn() if shuffle < 0.3 else rand())*shuffle)*dy*b, 
     827                      (iz+0.0+(randn() if shuffle < 0.3 else rand())*shuffle)*dz*c) 
     828               .rotate(*((randn(3) if rotate < 30 else rand(3))*rotate)) 
     829               for ix in range(nx) 
     830               for iy in range(ny) 
     831               for iz in range(nz)] 
     832    lattice = Composite(corners + faces_a + faces_b + faces_c) 
     833    return lattice 
    767834 
    768835SHAPE_FUNCTIONS = OrderedDict([ 
     
    775842]) 
    776843SHAPES = list(SHAPE_FUNCTIONS.keys()) 
     844LATTICE_FUNCTIONS = OrderedDict([ 
     845    ("sc", build_sc_lattice), 
     846    ("bcc", build_bcc_lattice), 
     847    ("fcc", build_fcc_lattice), 
     848]) 
     849LATTICE_TYPES = list(LATTICE_FUNCTIONS.keys()) 
    777850 
    778851def check_shape(title, shape, fn=None, show_points=False, 
     
    783856    r = shape.r_bins(q, r_step=r_step) 
    784857    sampling_density = samples / shape.volume 
     858    print("sampling points") 
    785859    rho, points = shape.sample(sampling_density) 
     860    print("calculating Pr") 
    786861    t0 = time.time() 
    787862    Pr = calc_Pr(r, rho-rho_solvent, points) 
     
    806881    Qx, Qy = np.meshgrid(qx, qy) 
    807882    sampling_density = samples / shape.volume 
     883    print("sampling points") 
    808884    t0 = time.time() 
    809885    rho, points = shape.sample(sampling_density) 
     
    844920                        help='lattice size') 
    845921    parser.add_argument('-z', '--spacing', type=str, default='2,2,2', 
    846                         help='lattice spacing') 
     922                        help='lattice spacing (relative to shape)') 
     923    parser.add_argument('-t', '--type', choices=LATTICE_TYPES, 
     924                        default=LATTICE_TYPES[0], 
     925                        help='lattice type') 
    847926    parser.add_argument('-r', '--rotate', type=float, default=0., 
    848927                        help="rotation relative to lattice, gaussian < 30 degrees, uniform otherwise") 
     
    860939    shuffle, rotate = opts.shuffle, opts.rotate 
    861940    shape, fn, fn_xy = SHAPE_FUNCTIONS[opts.shape](**pars) 
     941    view = tuple(float(v) for v in opts.view.split(',')) 
    862942    if nx > 1 or ny > 1 or nz > 1: 
    863         shape = build_cubic_lattice(shape, nx, ny, nz, dx, dy, dz, shuffle, rotate) 
     943        print("building %s lattice"%opts.type) 
     944        lattice = LATTICE_FUNCTIONS[opts.type] 
     945        shape = lattice(shape, nx, ny, nz, dx, dy, dz, shuffle, rotate) 
     946        # If comparing a sphere in a cubic lattice, compare against the 
     947        # corresponding paracrystalline model. 
     948        if opts.shape == "sphere" and dx == dy == dz: 
     949            radius = pars.get('radius', DEFAULT_SPHERE_RADIUS) 
     950            model_name = opts.type + "_paracrystal" 
     951            model_pars = { 
     952                "scale": 1., 
     953                "background": 0., 
     954                "lattice_spacing": 2*radius*dx, 
     955                "lattice_distortion": shuffle, 
     956                "radius": radius, 
     957                "sld": pars.get('rho', DEFAULT_SPHERE_CONTRAST), 
     958                "sld_solvent": 0., 
     959                "theta": view[0], 
     960                "phi": view[1], 
     961                "psi": view[2], 
     962            } 
     963            fn, fn_xy = wrap_sasmodel(model_name, **model_pars) 
    864964    title = "%s(%s)" % (opts.shape, " ".join(opts.pars)) 
    865965    if opts.dim == 1: 
     
    867967                    mesh=opts.mesh, qmax=opts.qmax, samples=opts.samples) 
    868968    else: 
    869         view = tuple(float(v) for v in opts.view.split(',')) 
    870969        check_shape_2d(title, shape, fn_xy, view=view, show_points=opts.plot, 
    871970                       mesh=opts.mesh, qmax=opts.qmax, samples=opts.samples) 
Note: See TracChangeset for help on using the changeset viewer.