Simulate X-ray Diffraction from Polycrystals in 3D.

The X-R ay D iffraction SIMULATOR package defines polycrystals as a mesh of tetrahedral single crystals and simulates diffraction as collected by a 2D discretized detector array while the sample is rocked around an arbitrary rotation axis. The full journal paper associated to the release of this code can be found here:

xrd_simulator: 3D X-ray diffraction simulation software supporting 3D polycrystalline microstructure morphology descriptions Henningsson, A. & Hall, S. A. (2023). J. Appl. Cryst. 56, 282-292. https://doi.org/10.1107/S1600576722011001

xrd_simulator was originally developed with the hope to answer questions about measurement optimization in scanning x-ray diffraction experiments. However, xrd_simulator can simulate a wide range of experimental diffraction setups. The essential idea is that the sample and beam topology can be arbitrarily specified, and their interaction simulated as the sample is rocked. This means that standard “non-powder” experiments such as scanning-3dxrd and full-field 3dxrd (or HEDM if you like) can be simulated as well as more advanced measurement sequences such as helical scans for instance. It is also possible to simulate powder like scenarios using orientation density functions as input.

Introduction

Before reading all the boring documentation (which is hosted here) let’s dive into some end to end examples to get us started on a good flavour.

The xrd_simulator is built around four python objects which reflect a diffraction experiment:

• A beam of x-rays (using the xrd_simulator.beam module)

• A 2D area detector (using the xrd_simulator.detector module)

• A 3D polycrystal sample (using the xrd_simulator.polycrystal module)

• A rigid body sample motion (using the xrd_simulator.motion module)

Once these objects are defined it is possible to let the detector collect scattering of the polycrystal as the sample undergoes the prescribed rigid body motion while being illuminated by the xray beam.

Let’s go ahead and build ourselves some x-rays:

import numpy as np
from xrd_simulator.beam import Beam

# The beam of xrays is represented as a convex polyhedron
# We specify the vertices in a numpy array.
beam_vertices = np.array(
   [
      [-1e6, -500.0, -500.0],
      [-1e6, 500.0, -500.0],
      [-1e6, 500.0, 500.0],
      [-1e6, -500.0, 500.0],
      [1e6, -500.0, -500.0],
      [1e6, 500.0, -500.0],
      [1e6, 500.0, 500.0],
      [1e6, -500.0, 500.0],
   ]
)

beam = Beam(
   beam_vertices,
   xray_propagation_direction=np.array([1.0, 0.0, 0.0]),
   wavelength=0.28523,
   polarization_vector=np.array([0.0, 1.0, 0.0]),
)

We will also need to define a detector:

from xrd_simulator.detector import Detector

# The detector plane is defined by it's corner coordinates det_corner_0,det_corner_1,det_corner_2
detector = Detector(
   det_corner_0=np.array([142938.3, -38400.0, -38400.0]),
   det_corner_1=np.array([142938.3, 38400.0, -38400.0]),
   det_corner_2=np.array([142938.3, -38400.0, 38400.0]),
   pixel_size=(55.0, 55.0),
   gaussian_sigma=1.0,
   max_gaussian_kernel_radius=5,
)

Next we go ahead and produce a sample, to do this we need to first define a mesh that describes the topology of the sample, in this example we make the sample shaped as a ball:

from xrd_simulator.mesh import TetraMesh

# xrd_simulator supports several ways to generate a mesh, here we
# generate meshed solid sphere using a level set.
mesh = TetraMesh.generate_mesh_from_levelset(
   level_set=lambda x: np.linalg.norm(x) - 768.0,
   bounding_radius=769.0,
   max_cell_circumradius=550.0,
)

Every element in the sample is composed of some material, or “phase”, we define the present phases in a list of xrd_simulator.phase.Phase objects, in this example only a single phase is present:

from xrd_simulator.phase import Phase

quartz = Phase(
   unit_cell=[4.926, 4.926, 5.4189, 90.0, 90.0, 120.0],
   sgname="P3221",  # (Quartz)
   path_to_cif_file=None,  # phases can be defined from crystalographic information files
)

The polycrystal sample can now be created. In this example the crystal elements have random orientations and the strain is uniformly zero in the sample:

from scipy.spatial.transform import Rotation as R
from xrd_simulator.polycrystal import Polycrystal

orientation = R.random(mesh.number_of_elements).as_matrix()
# element_phase_map assigns each mesh element to a phase (0 = first phase)
element_phase_map = np.zeros(mesh.number_of_elements, dtype=int)
polycrystal = Polycrystal(
   mesh,
   orientation,
   strain=np.zeros((3, 3)),
   phases=quartz,
   element_phase_map=element_phase_map,
)

We may save the polycrystal to disc by using the builtin save() command as

polycrystal.save('my_polycrystal', save_mesh_as_xdmf=True)

We can visualize the sample by loading the .xdmf file into your favorite 3D rendering program. In paraview the sampled colored by one of its Euler angles looks like this:

We can now define some motion of the sample over which to integrate the diffraction signal:

from xrd_simulator.motion import RigidBodyMotion

motion = RigidBodyMotion(
   rotation_axis=np.array([0, 1 / np.sqrt(2), -1 / np.sqrt(2)]),
   rotation_angle=np.radians(0.1),
   translation=np.array([123, -153.3, 3.42]),
)

Now that we have an experimental setup we may collect diffraction by letting the beam and detector interact with the sample:

peaks_dict = polycrystal.diffract(beam, motion, detector=detector)
diffraction_pattern, peaks_dict = detector.render(
   peaks_dict, frames_to_render=0, method="micro"
)

The resulting rendered detector frame will look something like the below. Note that the positions of the diffraction spots may vary as the crystal orientations were randomly generated!:

import matplotlib.pyplot as plt

fig, ax = plt.subplots(1, 1, figsize=(12, 12))
# render returns (frames, height, width), take first frame
pattern = (
   diffraction_pattern[0].cpu().numpy()
   if hasattr(diffraction_pattern, "cpu")
   else diffraction_pattern[0]
)
ax.imshow(pattern, cmap="gray", vmax=5000)
plt.show()

To compute several frames simply change the motion and collect the diffraction again. The sample may be moved before each computation using the same or another motion.

polycrystal.transform(motion, time=1.0)
peaks_dict = polycrystal.diffract(beam, motion, detector=detector)

Many more options for experimental setups and intensity rendering exist, have fun experimenting! The above example code can be found as a single .py file here.

Installation

Anaconda installation (Linux and Macos)

xrd_simulator is distributed on the conda-forge channel and the preferred way to install the xrd_simulator package is via Anaconda:

conda install -c conda-forge xrd_simulator
conda create -n xrd_simulator
conda activate xrd_simulator

This is meant to work across OS-systems and requires an Anaconda installation.

(The conda-forge feedstock of xrd_simulator can be found here.)

Anaconda installation (Windows)

xrd_simulator can be installed on Windows via Anaconda. The package now uses meshpy which provides better cross-platform support than the previous pygalmesh dependency.

Pip Installation

Pip installation is possible. The package uses meshpy which generally installs cleanly via pip:

pip install xrd-simulator

Source installation

Naturally one may also install from the sources:

git clone https://github.com/FABLE-3DXRD/xrd_simulator.git
cd xrd_simulator
python setup.py install

This will use meshpy which generally has better cross-platform support than pygalmesh.

Credits

xrd_simulator makes good use of xfab and meshpy for tetrahedral mesh generation. The source code of related repos can be found here:

• https://github.com/FABLE-3DXRD/xfab

• meshpy

Citation

If you feel that xrd_simulator was helpful in your research we would love for you to cite us.

xrd_simulator: 3D X-ray diffraction simulation software supporting 3D polycrystalline microstructure morphology descriptions Henningsson, A. & Hall, S. A. (2023). J. Appl. Cryst. 56, 282-292. https://doi.org/10.1107/S1600576722011001

Documentation

polycrystal

Polycrystal module for representing and simulating polycrystalline samples.

This module provides the Polycrystal class which handles:

• Multi-phase polycrystal representation as a tetrahedral mesh

• Diffraction computation

• Spatial transformations

• Crystal orientations and strains

class xrd_simulator.polycrystal.Polycrystal(mesh: TetraMesh, orientation: ndarray[Any, dtype[_ScalarType_co]] | Tensor, strain: ndarray[Any, dtype[_ScalarType_co]] | Tensor, phases: Phase | list[Phase], element_phase_map: ndarray[Any, dtype[_ScalarType_co]] | Tensor | None = None)

A multi-phase polycrystal represented as a tetrahedral mesh.

Each element in the mesh can be a single crystal. The polycrystal is created in laboratory coordinates with sample and lab coordinate systems initially aligned.

Parameters:

• mesh (TetraMesh) – Tetrahedral mesh defining sample geometry.

• orientation (numpy.ndarray or torch.Tensor) – Per-element orientation matrices (U) with shape (N, 3, 3) or (3, 3) if same for all elements.

• strain (numpy.ndarray or torch.Tensor) – Per-element Green-Lagrange strain tensor with shape (N, 3, 3) or (3, 3) if same for all elements.

• phases (Phase or list of Phase) – Single phase or list of phases present.

• element_phase_map (numpy.ndarray or torch.Tensor, optional) – Indices mapping elements to phases. Required for multi-phase samples. Default is None.

mesh_lab

Mesh in lab coordinates (updated during transforms).

Type:

TetraMesh

mesh_sample

Mesh in sample coordinates (fixed).

Type:

TetraMesh

orientation_lab

Crystal-to-lab orientation matrices.

Type:

torch.Tensor

orientation_sample

Crystal-to-sample orientation matrices.

Type:

torch.Tensor

strain_lab

Strain tensors in lab coordinates.

Type:

torch.Tensor

strain_sample

Strain tensors in sample coordinates.

Type:

torch.Tensor

phases

List of unique phases.

Type:

list of Phase

element_phase_map

Phase index per element.

Type:

torch.Tensor

diffract(beam: Beam, rigid_body_motion: RigidBodyMotion, min_bragg_angle: float = 0, max_bragg_angle: float = 1.5707963267948966, detector: Detector | None = None, verbose: bool = True) → Dict[str, Tensor | list]

Compute diffraction from the rotating and translating polycrystal.

Simulates diffraction while the sample is illuminated by an X-ray beam. Returns peaks that can be rendered on a detector.

Parameters:

• beam (Beam) – Monochromatic X-ray beam.

• rigid_body_motion (RigidBodyMotion) – Sample transformation over time domain [0, 1].

• min_bragg_angle (float, optional) – Minimum Bragg angle in radians. Default is 0.

• max_bragg_angle (float, optional) – Maximum Bragg angle in radians. Default is 90 * pi / 180.

• detector (Detector, optional) – Optional detector for automatic Bragg angle bounds.

• verbose (bool, optional) – If True, print progress messages during diffraction. Default is True.

Returns:

Dictionary containing:

• 'peaks': Tensor of diffraction peaks

• 'columns': Column names for peaks tensor

• 'mesh_lab': Mesh for computing convex hulls (if needed)

• 'rigid_body_motion': Motion object (if needed)

• 'beam': Beam object (if needed)

Return type:

dict

transform(rigid_body_motion, time)

Transform the polycrystal by performing a rigid body motion.

This function updates the polycrystal mesh (update in lab frame) with any dependent quantities such as face normals. Likewise, it updates the per-element crystal orientation matrices (U) as well as the lab frame description of strain tensors.

Parameters:

• rigid_body_motion (RigidBodyMotion) – Rigid body motion object describing the polycrystal transformation as a function of time on the domain time=[0, 1].

• time (float) – Time between [0, 1] at which to call the rigid body motion.

save(path: str, save_mesh_as_xdmf: bool = True) → None

Save polycrystal to disc via pickling.

Parameters:

• path (str) – File path at which to save, ending with the desired filename. The .pc extension is added if not present.

• save_mesh_as_xdmf (bool, optional) – If True, saves the polycrystal mesh with associated strains and crystal orientations as a .xdmf for visualization (sample coordinates). The results can be visualized with for instance ParaView (https://www.paraview.org/). The resulting data fields are the 6 unique components of the strain tensor (in sample coordinates) and the 3 Bunge Euler angles (Bunge, H. J. (1982). Texture Analysis in Materials Science. London: Butterworths.). Additionally, a single field specifying the material phases of the sample will be saved. Default is True.

classmethod load(path: str) → Polycrystal

Load polycrystal from disc via pickling.

Parameters:

path (str) – File path at which to load, ending with the desired filename.

Returns:

Loaded polycrystal object.

Return type:

Polycrystal

Raises:

ValueError – If the file does not end with .pc.

Warning

This function will unpickle data from the provided path. The pickle module is not intended to be secure against erroneous or maliciously constructed data. Never unpickle data received from an untrusted or unauthenticated source.

mesh

The mesh module is used to represent the morphology of a polycrystalline sample.

Once created and linked to a polycrystal, the mesh can be accessed directly through the xrd_simulator.polycrystal.Polycrystal.

Examples

Here is a minimal example of how to instantiate a mesh and save it to disc:

import numpy as np
import os
from xrd_simulator.mesh import TetraMesh

# Generate mesh with 2 elements from nodal coordinates.
nodal_coordinates = np.array([[0,0,0],[0,1,0],[1,0,0],[0,0,1],[0,0,-1]])
element_node_map  = np.array([[0,1,2,3],[0,1,2,4]])
mesh = TetraMesh.generate_mesh_from_vertices(nodal_coordinates, element_node_map)

# The mesh may be saved to disc for later usage or visualization.
artifacts_dir = os.path.join(os.path.dirname(__file__), 'test_artifacts')
os.makedirs(artifacts_dir, exist_ok=True)
mesh.save(os.path.join(artifacts_dir, 'my_mesh'))
mesh_loaded_from_disc = TetraMesh.load(os.path.join(artifacts_dir, 'my_mesh.xdmf'))

This should look something like this in a 3D viewer like paraview:

Below follows a detailed description of the mesh class attributes and functions.

class xrd_simulator.mesh.TetraMesh

Defines a 3D tetrahedral mesh with associated geometry data.

For level-set mesh generation, TetraMesh uses marching cubes from scikit-image to extract the surface, then meshpy (TetGen wrapper) to generate the volume tetrahedral mesh.

coord

Nodal coordinates, shape (n_nodes, 3). Each row defines the coordinates of a mesh node.

Type:

torch.Tensor

enod

Tetra element nodes, shape (n_elm, n_nodes). enod[i, :] gives the nodal indices of element i.

Type:

torch.Tensor

dof

Per node degrees of freedom. dof[i, :] gives the degrees of freedom of node i.

Type:

torch.Tensor

efaces

Element faces nodal indices, shape (n_elm, n_faces, 3). efaces[i, j, :] gives the nodal indices of face j of element i.

Type:

torch.Tensor

enormals

Element faces outward normals, shape (n_elm, n_faces, 3). enormals[i, j, :] gives the normal of face j of element i.

Type:

torch.Tensor

ecentroids

Per element centroids, shape (n_elm, 3).

Type:

torch.Tensor

eradius

Per element bounding ball radius, shape (n_elm,).

Type:

torch.Tensor

espherecentroids

Per element bounding ball centroids, shape (n_elm, 3).

Type:

torch.Tensor

evolumes

Per element volume, shape (n_elm,).

Type:

torch.Tensor

centroid

Global centroid of the entire mesh, shape (3,).

Type:

torch.Tensor

number_of_elements

Number of tetrahedral elements in the mesh.

Type:

int

classmethod generate_mesh_from_vertices(coord, enod)

Generate a mesh from vertices using the meshio package.

Parameters:

• coord (numpy.ndarray) – Nodal coordinates, shape (n_nodes, 3). Each row defines the coordinates of a mesh node.

• enod (numpy.ndarray) – Tetra element nodes, shape (n_elm, n_nodes). enod[i, :] gives the nodal indices of element i.

Returns:

The generated tetrahedral mesh.

Return type:

TetraMesh

classmethod generate_mesh_from_levelset(level_set, bounding_radius, max_cell_circumradius)

Generate a mesh from a level set using marching cubes + meshpy.

This method uses a two-step process:

1. Extract the surface from the level set using marching cubes (skimage)

2. Generate a volume tetrahedral mesh from the surface using meshpy

Parameters:

• level_set (callable) – Level set function. level_set(x) should give a negative output on the exterior of the mesh and positive on the interior.

• bounding_radius (float) – Bounding radius of mesh.

• max_cell_circumradius (float) – Upper bound for element radii.

Returns:

The generated tetrahedral mesh.

Return type:

TetraMesh

Raises:

ValueError – If marching cubes fails to extract a surface or if tetrahedral mesh generation fails.

translate(translation_vector)

Translate the mesh.

Parameters:

translation_vector (numpy.ndarray or torch.Tensor) – Translation vector [x, y, z], shape (3,).

rotate(rotation_axis, angle)

Rotate the mesh.

Parameters:

• rotation_axis (numpy.ndarray) – Rotation axis, shape (3,).

• angle (float) – Rotation angle in radians.

update(rigid_body_motion, time)

Apply a rigid body motion transformation to the mesh.

Parameters:

• rigid_body_motion (xrd_simulator.motion.RigidBodyMotion) – Rigid body motion object describing the polycrystal transformation as a function of time on the domain time=[0, 1].

• time (float) – Time between [0, 1] at which to call the rigid body motion.

save(file, element_data=None)

Save the tetra mesh to .xdmf paraview readable format.

Parameters:

• file (str) – Absolute path to save the mesh in .xdmf format.

• element_data (dict, optional) – Data associated to the elements. Keys are data names and values are lists or numpy arrays.

classmethod load(path)

Load a mesh from a saved mesh file using the meshio package.

Parameters:

path (str) – Absolute path to the mesh file.

Returns:

The loaded tetrahedral mesh.

Return type:

TetraMesh

phase

The phase module is used to represent material phase.

Each element of the xrd_simulator.polycrystal.Polycrystal mesh is linked to a xrd_simulator.phase.Phase describing things like the lattice reference unit cell and generating Miller indices of scattering planes. By providing a crystallographic information file (.cif), structure factors can be computed.

Examples

Here is a minimal example of how to instantiate a Phase object:

from xrd_simulator.phase import Phase

# The phase can be specified via a crystalographic information file.
quartz = Phase(unit_cell=[4.926, 4.926, 5.4189, 90., 90., 120.],
               sgname='P3221',
               path_to_cif_file='quartz.cif'
               )

The .cif file used in the above example can be found here.

Below follows a detailed description of the Phase class attributes and functions.

class xrd_simulator.phase.Phase(unit_cell, sgname, path_to_cif_file=None)

Defines properties related to a crystal class.

The barebone Phase object holds the space group name and unit cell of the crystal. From this it is possible to compute a set of Miller indices for a given wavelength that can give diffraction. In addition, it is possible to compute unit cell structure factors for the hkls which is useful to model the scattered intensity. This will however require a CIF file.

Parameters:

• unit_cell (list of float) – Crystal unit cell representation of the form [a, b, c, alpha, beta, gamma], where alpha, beta and gamma are in units of degrees while a, b and c are in units of angstrom.

• sgname (str) – Name of space group, e.g. 'P3221' for quartz (SiO2).

• path_to_cif_file (str, optional) – Path to CIF file. Default is None, in which case no structure factors are computed (i.e., structure_factors=None).

unit_cell

Crystal unit cell representation of the form [a, b, c, alpha, beta, gamma], where alpha, beta and gamma are in units of degrees while a, b and c are in units of angstrom.

Type:

list of float

sgname

Name of space group, e.g. 'P3221' for quartz (SiO2).

Type:

str

miller_indices

Allowable integer Miller indices (h, k, l) of shape (n, 3).

Type:

numpy.ndarray

structure_factors

Structure factors of allowable Miller indices (miller_indices) of shape (n, 2). structure_factors[i, 0] gives the real structure factor of hkl=miller_indices[i, :] while structure_factors[i, 1] gives the corresponding imaginary part.

Type:

numpy.ndarray

path_to_cif_file

Path to CIF file.

Type:

str

beam

The beam module is used to represent a beam of X-rays.

The idea is to create a xrd_simulator.beam.Beam object and pass it along to the xrd_simulator.polycrystal.Polycrystal.diffract() function to compute diffraction from the sample for the specified X-ray beam geometry.

Examples

Here is a minimal example of how to instantiate a beam object and save it to disc:

import numpy as np
import os
from xrd_simulator.beam import Beam

# The beam of xrays is represented as a convex polyhedron
# We specify the vertices in a numpy array.
beam_vertices = np.array([
    [-5., 0., 0.],
    [-5., 1., 0.],
    [-5., 0., 1.],
    [-5., 1., 1.],
    [5., 0., 0.],
    [5., 1., 0.],
    [5., 0., 1.],
    [5., 1., 1.]])

# The xray beam object is instantiated
beam = Beam(
    beam_vertices,
    xray_propagation_direction=np.array([1., 0., 0.]),
    wavelength=0.28523,
    polarization_vector=np.array([0., 1., 0.]))

# The xray beam may be saved to disc for later usage.
artifacts_dir = os.path.join(os.path.dirname(__file__), 'test_artifacts')
os.makedirs(artifacts_dir, exist_ok=True)
beam.save(os.path.join(artifacts_dir, 'my_xray_beam'))
beam_loaded_from_disc = Beam.load(os.path.join(artifacts_dir, 'my_xray_beam.beam'))

Below follows a detailed description of the beam class attributes and functions.

class xrd_simulator.beam.Beam(beam_vertices: Tensor | ndarray | list | tuple, xray_propagation_direction: Tensor | ndarray | list | tuple, wavelength: float, polarization_vector: Tensor | ndarray | list | tuple)

Represents a monochromatic X-ray beam as a convex polyhedra.

The beam has uniform intensity and is described in the laboratory coordinate system.

Parameters:

• beam_vertices (torch.Tensor or numpy.ndarray) – Vertices of the X-ray beam in units of microns, shape (N, 3).

• xray_propagation_direction (torch.Tensor or numpy.ndarray) – Propagation direction of X-rays, shape (3,).

• wavelength (float) – X-ray wavelength in units of angstrom.

• polarization_vector (torch.Tensor or numpy.ndarray) – Beam linear polarization unit vector, shape (3,). Must be orthogonal to the X-ray propagation direction.

vertices

Vertices of the X-ray beam in units of microns, shape (N, 3).

Type:

torch.Tensor

wavelength

X-ray wavelength in units of angstrom.

Type:

float

wave_vector

Beam wavevector with norm 2*pi/wavelength, shape (3,).

Type:

torch.Tensor

polarization_vector

Beam linear polarization unit vector, shape (3,). Must be orthogonal to the X-ray propagation direction.

Type:

torch.Tensor

centroid

Beam convex hull centroid, shape (3,).

Type:

torch.Tensor

halfspaces

Beam halfspace equation coefficients, shape (N, 4). A point x is on the interior of the halfspace if: halfspaces[i, :-1].dot(x) + halfspaces[i, -1] <= 0.

Type:

torch.Tensor

set_beam_vertices(beam_vertices: Tensor)

Set beam vertices and update all dependent quantities.

Parameters:

beam_vertices (torch.Tensor) – Vertices of the X-ray beam in units of microns, shape (N, 3).

save(path: str)

Save the X-ray beam to disc via pickling.

Parameters:

path (str) – File path at which to save, ending with the desired filename. The .beam extension is added if not present.

classmethod load(path: str) → Beam

Load the X-ray beam from disc via pickling.

Parameters:

path (str) – File path at which to load, ending with the desired filename.

Returns:

Loaded Beam object.

Return type:

Beam

Raises:

ValueError – If the file does not end with .beam.

Warning

This function will unpickle data from the provided path. The pickle module is not intended to be secure against erroneous or maliciously constructed data. Never unpickle data received from an untrusted or unauthenticated source.

detector

The detector module is used to represent a 2D area detector. After diffraction from a xrd_simulator.polycrystal.Polycrystal has been computed, the detector can render the scattering as a pixelated image via the xrd_simulator.detector.Detector.render() function.

Here is a minimal example of how to instantiate a detector object and save it to disc:

Examples:

import numpy as np
import os
from xrd_simulator.detector import Detector

# The detector is defined by it's corner coordinates det_corner_0, det_corner_1, det_corner_2
# Option 1: Specify number of pixels directly
detector = Detector(det_corner_0=np.array([142938.3, -38400., -38400.]),
                    det_corner_1=np.array([142938.3, 38400., -38400.]),
                    det_corner_2=np.array([142938.3, -38400., 38400.]),
                    n_pixels=(1024, 1396))  # (n_z, n_y) pixels

# Option 2: Specify pixel size - single value for square pixels
detector = Detector(det_corner_0=np.array([142938.3, -38400., -38400.]),
                    det_corner_1=np.array([142938.3, 38400., -38400.]),
                    det_corner_2=np.array([142938.3, -38400., 38400.]),
                    pixel_size=75.0)  # 75 µm square pixels

# Option 3: Specify pixel size - tuple for rectangular pixels
detector = Detector(det_corner_0=np.array([142938.3, -38400., -38400.]),
                    det_corner_1=np.array([142938.3, 38400., -38400.]),
                    det_corner_2=np.array([142938.3, -38400., 38400.]),
                    pixel_size=(75.0, 55.0))  # (z, y) in µm

# The detector may be saved to disc for later usage.
artifacts_dir = os.path.join(os.path.dirname(__file__), 'test_artifacts')
os.makedirs(artifacts_dir, exist_ok=True)
detector.save(os.path.join(artifacts_dir, 'my_detector'))
detector_loaded_from_disc = Detector.load(os.path.join(artifacts_dir, 'my_detector.det'))

Below follows a detailed description of the detector class attributes and functions.

class xrd_simulator.detector.Detector(det_corner_0: ndarray[Any, dtype[_ScalarType_co]], det_corner_1: ndarray[Any, dtype[_ScalarType_co]], det_corner_2: ndarray[Any, dtype[_ScalarType_co]], n_pixels: int | tuple | None = None, pixel_size: float | tuple | None = None, gaussian_sigma: float = 1.0, kernel_threshold: float = 0.02, max_gaussian_kernel_radius: int = 2, use_lorentz: bool = True, use_polarization: bool = True, use_structure_factor: bool = True)

Represents a rectangular 2D area detector.

The detector collects xrd_simulator.scattering_unit.ScatteringUnit during diffraction from a xrd_simulator.polycrystal.Polycrystal. The detector implements various rendering of the scattering as 2D pixelated frames. The detector geometry is described by specifying the locations of three detector corners. The detector is described in the laboratory coordinate system.

Parameters:

• det_corner_0 (numpy.ndarray) – Detector corner 3D coordinates, shape (3,). The origin of the detector is at this corner.

• det_corner_1 (numpy.ndarray) – Detector corner 3D coordinates, shape (3,). Defines the y-axis direction of the detector.

• det_corner_2 (numpy.ndarray) – Detector corner 3D coordinates, shape (3,). Defines the z-axis direction of the detector.

• n_pixels (int or tuple, optional) – Number of pixels. Can be a single int (square detector) or tuple (n_z, n_y). If provided, pixel sizes are calculated automatically from detector dimensions. Mutually exclusive with pixel_size.

• pixel_size (float or tuple, optional) – Pixel side length in microns. Can be a single value (square pixels) or tuple (pixel_size_z, pixel_size_y) for rectangular pixels. zdhat is the unit vector from det_corner_0 towards det_corner_2. ydhat is the unit vector from det_corner_0 towards det_corner_1.

• gaussian_sigma (float, optional) – Standard deviation of the Gaussian point spread function in pixels. Default is 1.0.

• kernel_threshold (float, optional) – Intensity threshold used to determine the kernel size. The kernel will extend until the Gaussian tail drops below this value. Default is 0.02.

• max_gaussian_kernel_radius (int) – Maximum radius of the Gaussian kernel. Default is 2 providing a 5x5 kernel, which is consistent with the Airy method. When the Gaussian kernel radius is larger than this value, the kernel will be truncated.

• use_lorentz (bool, optional) – Whether to apply the Lorentz factor correction. Default is True.

• use_polarization (bool, optional) – Whether to apply the polarization factor correction. Default is True.

• use_structure_factor (bool, optional) – Whether to apply the structure factor. Default is True.

pixel_size_z

Pixel side length along zdhat (rectangular pixels) in microns.

Type:

float

pixel_size_y

Pixel side length along ydhat (rectangular pixels) in microns.

Type:

float

det_corner_0

Detector corner 3D coordinates, shape (3,). The origin of the detector.

Type:

torch.Tensor

det_corner_1

Detector corner 3D coordinates, shape (3,).

Type:

torch.Tensor

det_corner_2

Detector corner 3D coordinates, shape (3,).

Type:

torch.Tensor

frames

Analytical diffraction frames.

Type:

list

zdhat

Detector basis vector along the z-axis direction.

Type:

torch.Tensor

ydhat

Detector basis vector along the y-axis direction.

Type:

torch.Tensor

normal

Detector normal, formed as the cross product: torch.cross(self.zdhat, self.ydhat).

Type:

torch.Tensor

zmax

Detector width along zdhat.

Type:

torch.Tensor

ymax

Detector height along ydhat.

Type:

torch.Tensor

pixel_coordinates

Real space 3D detector pixel coordinates, shape (n_z, n_y, 3).

Type:

torch.Tensor

gaussian_sigma

Standard deviation of the Gaussian point spread function.

Type:

float

kernel_threshold

Threshold value that determines the kernel size.

Type:

float

gaussian_kernel

Pre-computed normalized 2D Gaussian kernel.

Type:

torch.Tensor

lorentz_factor

Whether the Lorentz factor correction is applied.

Type:

bool

polarization_factor

Whether the polarization factor correction is applied.

Type:

bool

structure_factor

Whether the structure factor is applied.

Type:

bool

save(path: str) → None

Save detector to disk.

Parameters:

path (str) – Output file path. The .det extension is added if missing.

classmethod load(path: str) → Detector

Load detector from disk.

Parameters:

path (str) – Path to .det file.

Returns:

Loaded Detector instance.

Return type:

Detector

Raises:

ValueError – If file extension is not .det.

render(peaks_dict: Dict[str, Tensor | List[str]], frames_to_render: int = 1, method: str = 'auto', macro_grain_limit=None, micro_grain_limit=0.0010000000000000002, render_dtype: dtype | None = None, verbose: bool | None = None) → tuple[Tensor, Dict[str, Tensor | List[str]]]

Render diffraction frames using different peak profile methods.

Three physically motivated rendering methods are provided, each targeting a different grain-size regime. See the individual rendering methods for detailed physical justification.

• 'nano': Airy-disk profiles for sub-micron crystallites where the finite-lattice assumption breaks down and the Fourier transform of the lattice is no longer a delta. See _render_airy_peaks().

• 'micro': Gaussian point-spread for medium grains (~0.1–10 µm) where the instrument response dominates the spot shape. See _render_gauss_peaks().

• 'macro': Volume projection for grains larger than the pixel size, capturing crystal morphology and strain. See _render_projected_volumes().

• 'auto': Automatically selects per-peak based on grain volume (macro_grain_limit / micro_grain_limit thresholds).

Parameters:

• peaks_dict (dict) – Peak data dictionary containing 'peaks' tensor and metadata columns.

• frames_to_render (int, optional) – Number of frames to generate. The diffraction time domain [0, 1] is divided into frames_to_render equal bins, and each peak is assigned to the frame corresponding to its diffraction time. When set to 1, all peaks are rendered into a single integrated frame. Values less than 1 are clamped to 1. Default is 1.

• method (str, optional) – Rendering method: 'micro', 'nano', 'macro', or 'auto'. Default is 'micro'.

• macro_grain_limit (float, optional) – Volume threshold (cubic microns) above which grains use the macro rendering method. Default is pixel_size**3.

• micro_grain_limit (float, optional) – Volume threshold (cubic microns) below which grains use the nano rendering method. Default is 0.001 (0.1³ cubic microns).

• render_dtype (torch.dtype, optional) – Data type for rendered frames. If None, uses the dtype of the peaks tensor.

• verbose (bool or None, optional) – If True, print progress messages during rendering. If None, inherits the verbose flag from peaks_dict (set by polycrystal.diffract). Defaults to True if neither is set.

Returns:

• diffraction_frames (torch.Tensor) – Rendered diffraction frames with shape (frames, height, width).

• peaks_dict (dict) – Updated peak data dictionary with detector intersection columns (zd, yd, incident_angle, frame, intensity_factors) appended and non-visible / invalid peaks filtered out. Identical to the output of peaks_detector_intersection().

Raises:

ValueError – If an invalid method is specified.

contains(zd: Tensor, yd: Tensor) → Tensor

Check if detector coordinates are within bounds.

Parameters:

• zd (torch.Tensor) – Z-axis detector coordinates.

• yd (torch.Tensor) – Y-axis detector coordinates.

Returns:

Boolean tensor indicating which coordinates are within detector bounds.

Return type:

torch.Tensor

peaks_detector_intersection(peaks_dict: Dict[str, Tensor | List[str]], frames_to_render: int = 1, verbose: bool = True) → Dict[str, Tensor | List[str]]

Project diffraction peaks onto the detector and filter to those that hit it.

Takes the peaks_dict produced by Polycrystal.diffract() and determines where each diffracted ray intersects the detector plane. Peaks that fall outside the detector area or that have physically invalid intensity (e.g. infinite Lorentz factor) are discarded.

The method performs three steps in order:

1. Intersection — computes where each diffracted ray (defined by K_out and Source columns) hits the detector plane, appending zd, yd and incident_angle.

2. Filtering — discards peaks that miss the detector (via contains()) or have non-finite / non-positive intensity factors.

3. Frame assignment & intensity — bins peaks into frames_to_render equal time bins and computes the combined intensity_factors column from the enabled correction factors (structure, polarization, Lorentz).

The following columns are appended to the peaks tensor:

Col	Name	Description
25	zd	Detector z-coordinate of the hit (µm).
26	yd	Detector y-coordinate of the hit (µm).
27	incident_angle	Angle between the diffracted ray and the detector normal (degrees).
28	frame	Frame index assigned from the diffraction time via equal-width binning over [0, 1].
29	intensity_factors	Product of the enabled correction factors (structure, polarization, Lorentz) giving the scattering strength per unit volume.

Parameters:

• peaks_dict (dict) – Peak data dictionary as returned by Polycrystal.diffract(). Must contain at least 'peaks' (torch.Tensor) and 'columns' (list of str).

• frames_to_render (int, optional) – Number of time frames. The diffraction-time interval [0, 1] is split into frames_to_render equal bins and each peak is assigned to the bin its diffraction_time falls into. Default is 1.

• verbose (bool, optional) – If True, emit warnings when peaks with invalid intensity factors are discarded. Default is True.

Returns:

A new dictionary with the same keys as peaks_dict but with 'peaks' and 'columns' updated (augmented columns, filtered rows). The original peaks_dict is not modified.

Return type:

dict

Examples

>>> peaks_dict = polycrystal.diffract(beam, motion)
>>> result = detector.peaks_detector_intersection(peaks_dict,
...                                              frames_to_render=5)
>>> result["columns"][-5:]
['zd', 'yd', 'incident_angle', 'frame', 'intensity_factors']

motion

The motion module is used to represent a rigid body motion.

During diffraction from a xrd_simulator.polycrystal.Polycrystal, the xrd_simulator.motion.RigidBodyMotion object describes how the sample is translating and rotating. The motion can be used to update the polycrystal position via the xrd_simulator.polycrystal.Polycrystal.transform() function.

Examples

Here is a minimal example of how to instantiate a rigid body motion object, apply the motion to a pointcloud and save the motion to disc:

import numpy as np
import os
from xrd_simulator.motion import RigidBodyMotion

# The motion is described by a rotation axis and a translation vector
motion = RigidBodyMotion(rotation_axis=np.array([0, 1/np.sqrt(2), -1/np.sqrt(2)]),
                         rotation_angle=np.radians(2.0),
                         translation=np.array([123, -153.3, 3.42]))

# The motion can be applied to transform a set of points
# Points should have shape (N, 3) where N is the number of points
points = np.random.rand(10, 3)
transformed_points = motion( points, time=0.421 )

# The motion may be saved to disc for later usage.
artifacts_dir = os.path.join(os.path.dirname(__file__), 'test_artifacts')
os.makedirs(artifacts_dir, exist_ok=True)
motion.save(os.path.join(artifacts_dir, 'my_motion'))
motion_loaded_from_disc = RigidBodyMotion.load(os.path.join(artifacts_dir, 'my_motion.motion'))

Below follows a detailed description of the RigidBodyMotion class attributes and functions.

class xrd_simulator.motion.RigidBodyMotion(rotation_axis, rotation_angle, translation, origin=tensor([0., 0., 0.]))

Rigid body transformation by Euler axis rotation and translation.

A rigid body motion is defined in the laboratory coordinates system.

The motion is parametric in the interval time=[0, 1] and will perform a rigid body transformation of a point x by linearly uniformly rotating it from [0, rotation_angle] and translating [0, translation]. I.e., if called at a time time=t, the motion will first rotate the point t * rotation_angle radians around rotation_axis and next translate the point by the vector t * translation.

Parameters:

• rotation_axis (numpy.ndarray or torch.Tensor) – Rotation axis, shape (3,).

• rotation_angle (float) – Radians for final rotation when time=1.

• translation (numpy.ndarray or torch.Tensor) – Translation vector, shape (3,).

• origin (numpy.ndarray or torch.Tensor, optional) – Point in space about which the rigid body motion is defined. Defaults to the origin (0, 0, 0). All translations are executed in relation to the origin and all rotations are rotations about the point of origin. Shape (3,).

rotation_axis

Rotation axis, shape (3,).

Type:

torch.Tensor

rotation_angle

Radians for final rotation when time=1.

Type:

torch.Tensor

translation

Translation vector, shape (3,).

Type:

torch.Tensor

origin

Point in space about which the rigid body motion is defined, shape (3,).

Type:

torch.Tensor

__call__(vectors, time)

Find the transformation of a set of points at a prescribed time.

Calling this method executes the rigid body motion with respect to the currently set origin.

Parameters:

• vectors (numpy.ndarray or torch.Tensor) – A set of points to be transformed, shape (3, N), (N, 3), or (N, 4, 3).

• time (float, numpy.ndarray, or torch.Tensor) – Time to compute for. Can be scalar or shape (N,) for per-vector times.

Returns:

Transformed vectors with shape (3, N), (N, 3), or (N, 4, 3).

Return type:

torch.Tensor

rotate(vectors, time)

Find the rotational transformation of a set of vectors.

This function only applies the rigid body rotation and will not respect the origin of the motion. This function is intended for rotation of diffraction and wavevectors. Use the __call__ method to perform a physical rigid body motion respecting the origin.

Parameters:

• vectors (numpy.ndarray or torch.Tensor) – A set of points in 3D Euclidean space to be rotated, shape (3, N) or (N, 3).

• time (float, numpy.ndarray, or torch.Tensor) – Time to compute for. Can be scalar or shape (N,) for per-vector times.

Returns:

Transformed vectors of shape (3, N) or (N, 3).

Return type:

torch.Tensor

translate(vectors, time)

Find the translational transformation of a set of points.

This function only applies the rigid body translation.

Parameters:

• vectors (numpy.ndarray or torch.Tensor) – A set of points in 3D Euclidean space to be translated, shape (3, N) or (N, 3).

• time (float) – Time to compute for.

Returns:

Transformed vectors of shape (3, N) or (N, 3).

Return type:

torch.Tensor

inverse()

Create an instance of the inverse motion.

The inverse is defined by negative translation and rotation axis vectors.

Returns:

The inverse motion with a reversed rotation and translation.

Return type:

RigidBodyMotion

save(path)

Save the motion object to disc via pickling.

Parameters:

path (str) – File path at which to save, ending with the desired filename. The .motion extension is added if not present.

classmethod load(path)

Load the motion object from disc via pickling.

Parameters:

path (str) – File path at which to load, ending with the desired filename.

Returns:

Loaded motion object.

Return type:

RigidBodyMotion

Raises:

ValueError – If the file does not end with .motion.

Warning

This function will unpickle data from the provided path. The pickle module is not intended to be secure against erroneous or maliciously constructed data. Never unpickle data received from an untrusted or unauthenticated source.

templates

Templates module for fast creation of sample types and diffraction geometries.

The templates module allows for fast creation of a few select sample types and diffraction geometries without having to worry about any of the “under the hood” scripting.

xrd_simulator.templates.s3dxrd(parameters)

Construct a scanning-three-dimensional X-ray diffraction experiment.

This is a helper/utility function for quickly creating an experiment. For full control over the diffraction geometry consider custom creation of the primitive quantities: xrd_simulator.beam.Beam, and xrd_simulator.detector.Detector separately.

Parameters:

parameters (dict) –

Dictionary with the following fields:

• "detector_distance": Distance from sample origin to detector centre in units of microns.

• "number_of_detector_pixels_z": Number of detector pixels along z-axis.

• "number_of_detector_pixels_y": Number of detector pixels along y-axis.

• "detector_center_pixel_z": Intersection pixel coordinate between beam centroid line and detector along z-axis.

• "detector_center_pixel_y": Intersection pixel coordinate between beam centroid line and detector along y-axis.

• "pixel_side_length_z": Detector pixel side length in units of microns along z-axis.

• "pixel_side_length_y": Detector pixel side length in units of microns along y-axis.

• "wavelength": Wavelength in units of Angstrom.

• "beam_side_length_z": Beam side length in units of microns.

• "beam_side_length_y": Beam side length in units of microns.

• "rotation_step": Angular frame integration step in units of radians.

• "rotation_axis": Axis around which to positively rotate the sample by rotation_step radians.

Returns:

Objects defining an experiment: (xrd_simulator.beam.Beam, xrd_simulator.detector.Detector, xrd_simulator.motion.RigidBodyMotion).

Return type:

tuple

Examples

import numpy as np
from xrd_simulator import templates

parameters = {
    "detector_distance": 191023.9164,
    "detector_center_pixel_z": 256.2345,
    "detector_center_pixel_y": 255.1129,
    "pixel_side_length_z": 181.4234,
    "pixel_side_length_y": 180.2343,
    "number_of_detector_pixels_z": 512,
    "number_of_detector_pixels_y": 512,
    "wavelength": 0.285227,
    "beam_side_length_z": 512 * 200.,
    "beam_side_length_y": 512 * 200.,
    "rotation_step": 1.0,
    "rotation_axis": np.array([0., 0., 1.0])
}

beam, detector, motion = templates.s3dxrd(parameters)

xrd_simulator.templates.polycrystal_from_odf(orientation_density_function, number_of_crystals, sample_bounding_cylinder_height, sample_bounding_cylinder_radius, unit_cell, sgname, path_to_cif_file=None, maximum_sampling_bin_seperation=0.08726646259971647, strain_tensor=<function <lambda>>)

Fill a cylinder with crystals from a given orientation density function.

The orientation_density_function is sampled by discretizing orientation space over the unit quaternions. Each bin is assigned its appropriate probability, assuming the orientation_density_function is approximately constant over a single bin. Each sampled orientation is constructed by first drawing a random bin and next drawing uniformly from within that bin, again assuming that orientation_density_function is approximately constant over a bin.

Parameters:

• orientation_density_function (callable) – Function orientation_density_function(x, q) -> float where input variable x is a numpy array of shape (3,) representing a spatial coordinate in the cylinder (x, y, z) and q is a numpy array of shape (4,) representing an orientation in SO3 by a unit quaternion. The format of the quaternion is “scalar last” (same as in scipy.spatial.transform.Rotation).

• number_of_crystals (int) – Approximate number of crystal elements to compose the cylinder volume.

• sample_bounding_cylinder_height (float) – Height of sample cylinder in units of microns.

• sample_bounding_cylinder_radius (float) – Radius of sample cylinder in units of microns.

• unit_cell (list of float) – Crystal unit cell representation of the form [a, b, c, alpha, beta, gamma], where alpha, beta and gamma are in units of degrees while a, b and c are in units of Angstrom.

• sgname (str) – Name of space group, e.g. 'P3221' for quartz, SiO2, for instance.

• path_to_cif_file (str, optional) – Path to CIF file. Default is None, in which case no structure factors are computed.

• maximum_sampling_bin_seperation (float, optional) – Discretization steplength of orientation space using spherical coordinates over the unit quaternions in units of radians. A smaller steplength gives more accurate sampling of the input orientation_density_function but is computationally slower. Default is np.radians(5.0).

• strain_tensor (callable, optional) – Strain tensor field over sample cylinder. strain_tensor(x) -> numpy.ndarray of shape (3, 3) where input variable x is a numpy array of shape (3,) representing a spatial coordinate in the cylinder (x, y, z). Default returns zero strain.

Returns:

A polycrystal sample filling the specified cylinder.

Return type:

Polycrystal

Examples

import numpy as np
from xrd_simulator.templates import polycrystal_from_odf


# uniform orientation distribution function.
def ODF(x, q): return 1. / (np.pi**2)


number_of_crystals = 500
bounding_height = 50.0
bounding_radius = 25.0
unit_cell = [4.926, 4.926, 5.4189, 90., 90., 120.]
sgname = 'P3221'  # Quartz
max_bin = np.radians(10.0)
path_to_cif_file = None

def strain_tensor(x): return np.array([[0, 0, 0.02 * x[2] / bounding_height],
                                       [0, 0, 0],
                                       [0, 0, 0]])  # Linear strain gradient along rotation axis.


polycrystal = polycrystal_from_odf(ODF,
                                   number_of_crystals,
                                   bounding_height,
                                   bounding_radius,
                                   unit_cell,
                                   sgname,
                                   path_to_cif_file,
                                   max_bin,
                                   strain_tensor)

More complicated ODFs are also possible to use. Why not use a von-Mises-Fisher distribution for instance:

from scipy.special import iv
import numpy as np
# Von Mises-Fisher orientation distribution function.
# https://en.wikipedia.org/wiki/Von_Mises%E2%80%93Fisher_distribution
def ODF(x, q):
    kappa = 0.064 * x.dot(x) + 0.001
    mu    = np.array([ 0.,  0.95689793,  0.28706938, -0.0440173 ])
    p     = 4.0
    I     = iv( p/2 - 1, kappa)
    Cp    = kappa**(p/2 - 1) / ( (2*np.pi)**(p/2) * I )
    return Cp * ( np.exp( kappa * np.dot( q, mu ) ) + np.exp( -(kappa * np.dot( q, mu )) ) )

xrd_simulator.templates.get_uniform_powder_sample(sample_bounding_radius, number_of_grains, unit_cell, sgname, strain_tensor=array([[0., 0., 0.], [0., 0., 0.], [0., 0., 0.]]), path_to_cif_file=None)

Generate a polycrystal with grains overlayed at origin, uniform orientations.

Parameters:

• sample_bounding_radius (float) – Bounding radius of sample. All tetrahedral crystal elements will be overlayed within a sphere of sample_bounding_radius radius.

• number_of_grains (int) – Number of grains composing the polycrystal sample.

• unit_cell (list of float) – Crystal unit cell representation of the form [a, b, c, alpha, beta, gamma], where alpha, beta and gamma are in units of degrees while a, b and c are in units of Angstrom.

• sgname (str) – Name of space group, e.g. 'P3221' for quartz, SiO2, for instance.

• strain_tensor (numpy.ndarray, optional) – Strain tensor to apply to all tetrahedral crystal elements contained within the sample, shape (3, 3). Default is zero strain.

• path_to_cif_file (str, optional) – Path to CIF file. Default is None, in which case no structure factors are computed.

Returns:

A polycrystal sample with number_of_grains grains.

Return type:

Polycrystal

Examples

"""Example: Create a uniform powder sample using templates.get_uniform_powder_sample.

This example demonstrates how to create a polycrystal sample with uniformly
distributed grain orientations, suitable for powder diffraction simulations.
"""
import numpy as np
from xrd_simulator import templates

# Define sample parameters
sample_bounding_radius = 100.0  # microns
number_of_grains = 500  # Number of randomly oriented grains

# Define phase parameters (quartz example)
unit_cell = [4.926, 4.926, 5.4189, 90., 90., 120.]  # a, b, c, alpha, beta, gamma
sgname = 'P3221'  # Space group name

# Optional: strain tensor (default is zero strain)
strain_tensor = np.zeros((3, 3))

# Create the polycrystal sample
polycrystal = templates.get_uniform_powder_sample(
    sample_bounding_radius=sample_bounding_radius,
    number_of_grains=number_of_grains,
    unit_cell=unit_cell,
    sgname=sgname,
    strain_tensor=strain_tensor,
    path_to_cif_file=None,  # Optional: path to CIF file for structure factors
)

print(f"Created polycrystal with {polycrystal.mesh_sample.number_of_elements} elements")
print(f"Number of phases: {len(polycrystal.phases)}")

laue

Collection of functions for solving time-dependent Laue equations.

This module provides functions for arbitrary rigid body motions. It is mainly used internally by the xrd_simulator.polycrystal.Polycrystal. However, for the advanced user, access to these functions may be of interest.