nanoBragg

program for calculation of absolute scattering from molecules and small crystals

This short program calculates the absolute-scale scattering from a nanocrystal that is "bathed" in a beam of a given integrated photon density (specified in photons/meter2). For example, 1012 photons focused into a 3-micron round beam is represented by "-fluence 1.4e24". Images of the expected photons/pixel on the detector, with and without photon-counting noise are generated in SMV format (suitable for display with adxv, MOSFLM, or most any other diffraction image display program).
The structure factor of the spots should be provided on an absolute "electron" scale (as output by programs like phenix.fmodel, REFMAC or SFALL), but must be converted to a plain text file of h,k,l,F. Note that no symmetry is imposed by this program, not even Friedel symmetry, so all reflections you wish to be non-zero intensity must be specified, including F000. The unit cell and crystal orientation may be provided as a MOSFLM-style orientation matrix, which is again a text file and the first nine tokens read from it are taken as the x,y,z components of the three reciprocal-space cell vectors (a,b,c are the columns, x,y,z are the rows).
The program also contains an option for adding approximate scattering from the water droplet presumed to be surrounding the nanocrystal. The diameter of this droplet in microns is provided with the "-water" option, and assumes a forward-scattering structure factor of 2.57 electrons. The default value for this option is zero.
source code: nanoBragg.c (49k) or binaries for linux, OSX, Cygwin, Windows CUDA
auxillary programs:
UBtoA.awk can be used to generate a MOSFLM -style orientation matrix, and
mtz_to_P1hkl.com is a script for converting mtz-formatted structure factors into a format that nanoBragg can read.
noisify is a program that takes the "photons/pixel" noiseless intensity values output by nonBragg or nanoBragg, or nearBragg as "floagimage.bin" and adds different kinds of noise to it to generate an SMV file. This is usually faster than re-running nonBragg just to change things like beam intensity. In addition to photon shot noise, noisify has a few kinds of noise that nonBragg doesn't implement, such as pixel read-out noise, beam flicker, and calibration error.
float_add may be used to add the raw "float" binary files output by nonBragg, nanoBragg, or even nearBragg so that renderings may be divided up on separate CPUs and then combined together. The resulting raw files may then be converted to SMV images with noisify.
float_func can perform a large number of operations on these "floagimage.bin" files.
nonBragg is for generating scattering from amorphous substances, like water and air. You will need to feed it a text file containing the "structure factor" of the amorphous material vs sin(theta)/lambda. A few examples are:
air.stol He.stol ice.stol nanoice.stol Paratone-N.stol water.stol

example usage:

compile it 

gcc -O -O -o nanoBragg nanoBragg.c -lm -static
get some structure factor data 
getcif.com 3pcq
refine to get Fs on an absolute scale 
refmac5 hklin 3pcq.mtz xyzin 3pcq.pdb hklout refmacout.mtz xyzout refmacout.pdb << EOF | tee refmac.log
REFI TYPE RIGID
TWIN
EOF
extract the (de-twinned) calculated Fs, which are always 100% complete: 
mtz_to_P1hkl.com refmacout.mtz FC_ALL_LS
make a random orientation matrix: 
./UBtoA.awk << EOF | tee A.mat
CELL 281 281 165.2 90 90 120 
WAVE 6.2
RANDOM
EOF
run the simulation of a 10x10x10 unit cell crystal 
./nanoBragg -hkl P1.hkl -matrix A.mat -lambda 6.2 -N 10
view the result 
adxv intimage.img
convert and re-scale as regular graphics file 
convert -depth 16 -type Grayscale -colorspace GRAY -endian LSB -size 1024x1024+512 \
  -negate -normalize \
  GRAY:intimage.img intimage.png
Note the "low resolution hole", which is due to the missing low-angle data in the PDB deposition. Missing high-resolution spots that would otherwise fall on the detector will generate a WARNING message in the program output and potentially undefined spot intensities, so make sure you "fill" the resolution of interest in the P1.hkl file.
Note also that this image is very clear, with lots of inter-Bragg spot subsidiary peaks. That is because it is a noiseless simulation. Now have a look at the "noiseimage.img" which is scaled so that one pixel unit is one photon:.
adxv noiseimage.img
It has been re-scaled here as a png for better viewing:

Still not bad, but this is because there is no background, and the default fluence is 1e24, or 1012 photons focused into a 1 micron beam.
Now lets do something more realistic. The fluence of a 1012-photon pulse focused into a 7 micron beam is 2e22 photons/m2. Also, the liquid jet used by Chapman et al (2010) was four microns wide: 
./nanoBragg -hkl P1.hkl -matrix A.mat -lambda 6.2 -N 10 -fluence 2e22 -water 4

adxv noiseimage.img

If you look closely, you can see the spots. Note that this is an idealized case where only photon-counting noise is present. there is no detector read-out noise, no point-spread function, no amplifier drift, no pixel saturation and no calibration errors. Many of these errors can be added using noisify.c, but not all. Watch this space for updates.

SAXS simulations

nanoBragg can also be used to simulate small-angle X-ray scattering (SAXS) patterns by simply setting the number of unit cells to one (-N 1 on the command line). Tricubic interpolation between the hkl indicies will be used to determine the intensity between the "spots".

example usage:

get lysozyme 

getcif.com 193l
refine to get the solvent parameters 
phenix.refine 193l.pdb 193l.mtz | tee phenix_refine.log
Now put these atoms into a very big unit cell. It is important that this cell be at least 3-4 times bigger than your molecule in all directions. Otherwise, you will get neigbor-interference effects. Most people don't want that in their SAXS patterns. 
pdbset xyzin 193l.pdb xyzout bigcell.pdb << EOF
CELL 250 250 250 90 90 90
SPACEGROUP 1
EOF
calculate structure factors of the molecule isolated in a huge "bath" of the best-fit solvent. 
phenix.fmodel bigcell.pdb high_resolution=10 \
 k_sol=0.35 b_sol=46.5 mask.solvent_radius=0.5 mask.shrink_truncation_radius=0.16
note that this procedure will fill the large cell with a solvent of average electron density 0.35 electrons/A^3. The old crystallographic contacts will be replaced with the same solvent boundary model that fit the solvent channels in the crystal structure.
now we need to convert these Fs into a format nanoBragg can read 
mtz_to_P1hkl.com bigcell.pdb.mtz
and create a random orientation matrix 
./UBtoA.awk << EOF | tee bigcell.mat
CELL 250 250 250 90 90 90 
WAVE 1
RANDOM
EOF
and now, make the diffraction image 
./nanoBragg -mat bigcell.mat -hkl P1.hkl -lambda 1 -dispersion 0 \
  -distance 1000 -detsize 100 -pixel 0.1 \
  -hdiv 0 -vdiv 0 \
  -fluence 1e32 -N 1
adxv noiseimage.img


Notice that the center of the image is white. This is not a beamstop! What is actually going on is that F000 is missing in P1.hkl, and so is being replaced with zero intensity. You can fix this by adding an F000 term: 
echo "0 0 0 520" >> P1.hkl
./nanoBragg -mat bigcell.mat -hkl P1.hkl -lambda 1 -dispersion 0 \
  -distance 1000 -detsize 100 -pixel 0.1 \
  -hdiv 0 -vdiv 0 \
  -fluence 1e32 -N 1
adxv noiseimage.img


You might wonder, however, why the F000 term is ~500 and not the number of electrons in lysozyme, which is ~8000. The reason here is the bulk solvent. The volume of water displaced by the lysozyme molecule contains almost as many electrons as the lysozyme molecule itself. Protein, however, is slightly denser, and so there are an extra ~520 electrons "peeking" above the average density of the solvent.
Of course, most real SAXS patterns are centrosymmetric because they are an average over trillions of molecules in solution, each in a random orientation. The SAXS pattern generated here is for a single molecule exposed to 1e34 photons/m2, but this is equivalent to 1e18 photons focused onto an area barely larger than the molecule! However, 1e12 photons focused onto a 100x100 micron area (1e20 phm2) containing 1e14 molecules will generate a pattern of similar intensity level, albeit rotationally averaged.
One way to simulate such images would be to use this program to generate a few thousand or million orientations and then average the results. This could be instructive for exploring fluctuation SAXS.

However, a much faster way would be to pre-average the squared structure factors to form a new P1.hkl file, and then generate one image with a "-fluence" equal to the actual fluence, multiplied by the number of exposed molecules. A convenient script for doing this is: 

mtz_to_stol.com bigcell.pdb.mtz
which will create a file called mtz.stol that you can feed to nonBragg: 
./nonBragg -stol mtz.stol -lambda 1 -dispersion 0 \
  -distance 1000 -detsize 100 -pixel 0.1 \
  -hdiv 0 -vdiv 0 \
  -flux 1e13 -thick 1.2 -MW 14000 -density 0.01
adxv noiseimage.img


Which starts to look more like a SAXS pattern from a conventional SAXS beamline. Note that the "density" of the sample in this case is 0.01 g/cm^3 or 10 mg/mL.
Command-line options:
-hkl filename.hkl
the structure factor text list. Default: re-read dumpfile from last run
-matrix auto.mat
cell/orientation matrix file, takes first nine text numbers found
-cell a b c alpha beta gamma
specify unit cell dimensions (Angstrom and degrees)
-misset
instead of matrix, specify MOSFLM-style misseting angles about x,y,z (degrees)
-Na
number of unit cells along crystal a axis
-Nb
number of unit cells along crystal b axis
-Nc
number of unit cells along crystal c axis
-N
number of unit cells in all three directions (ovverides above)
-samplesize
alternative: linear dimension of the crystal in all three directions (mm)
-sample_thick or -sample_x ; -sample_width or -sample_y ; -sample_height or -sample_z
alternative: linear dimension of the crystal in specified directions (mm)
-img filename.img
optional: inherit and interpret header of an existing SMV-format diffraction image file
-distance
distance from sample to beam center on detector (mm)
-close_distance
distance from sample to nearest point in detector plane, XDS-style (mm)
-detsize
detector size in x and y (mm)
-detsize_x
detector size in x direction (mm)
-detsize_y
detector size in y direction (mm)
-pixel
detector pixel size (mm)
-detpixels
detector size in x and y (pixels)
-detpixels_x
detector size in x direction (pixels)
-detpixels_y
detector size in y direction (pixels)
-Xbeam
direct beam position in x direction (mm) Default: center
-Ybeam
direct beam position in y direction (mm) Default: center
-Xclose -Yclose
instead of beam center, specify point on detector closest to the sample (mm) Default: derive from Xbeam Ybeam
-ORGX -ORGY
instead of beam center, specify XDS-stype point on detector closest to the sample (pixels) Default: derive from Xbeam Ybeam
-detector_rotx -detector_roty -detector_rotz
specify detector mis-orientation rotations about x,y,z axes (degrees)
-twotheta
specify detector rotation about sample (degrees)
-pivot sample|beam
specify if detector rotations should be about the crystal or about the beam center point on the detector surface
-xdet_vector -ydet_vector -zdet_vector -beam_vector -polar_vector -spindle_axis -twotheta_axis
explicity define unit vectors defining detector and beam orientation (XDS style)
-pix0_vector
explicity define XYZ coordinate of the first pixel in the output file (as printed in the output)
-curved_det
all detector pixels same distance from sample (origin)
-oversample
number of sub-pixels per pixel. Default: 1
-roi xmin xmax ymin ymax
only render pixels within a set range. Default: all detector
-mask mask.img
optional: skip over pixels that have zero value in a provided SMV-format image file
-lambda
incident x-ray wavelength (Angstrom). Default: 1
-fluence
incident x-ray intensity (photons/m^2). Default: 1.26e29 so I=F^2
-flux
incident x-ray intensity (photons/s). Default: none
-exposure
exposure time (s) used to convert flux and beam size to fluence. Default: 1
-beamsize
linear size of incident x-ray beam at sample (mm). Default: 0.1
-hdivrange
horizontal angular spread of source points (mrad). Default: 0
-vdivrange
vertical angular spread of source points (mrad). Default: 0
-hdivstep
number of source points in the horizontal. Default: 1
-vdivstep
number of source points in the vertical. Default: 1
-round_div -square_div
make the 2D divergence distribution round or square. Default: round
-dispersion
spectral dispersion: delta-lambda/lambda (percent). Default: 0
-dispsteps
number of wavelengths in above range. Default: 1
-sourcefile
optionally specify a text file containing x,y,z,relative_intensity,wavelength of each desired point source
-coherent
coherently add everything, even different wavelengths. Not the default
-mosaic
simulate mosaic spread with random points on a spherical cap of specified diameter (degrees). Default: 0
-mosaic_domains
number of discrete mosaic domains to render. Default: 10 if mosaic>0 recommend a lot more
-phi -osc -phistep or -phisteps
simulate a spindle rotation about the spindle axis by averaging a series of stills. Default: 0
-phistep
angular step for simulating phi spindle rotation (deg). Default: derive from phisteps
-phisteps
number of steps for simulating phi spindle rotation (). Default: 2 if osc>0 recommend a lot more
-floatfile
name of binary pixel intensity output file (4-byte floats)
-intfile
name of smv-formatted output file.
-pgmfile
name of pgm-formatted output file.
-noisefile
name of smv-formatted output file containing photon-counting noise.
-nonoise
do not calculate noise or output noisefile
-nopgm
do not output pgm file
-scale
scale factor for intfile. Default: fill dynamic range
-pgmscale
scale factor for the pgm output file. Default: fill dynamic range
-adcoffset
specify the zero-photon level in the output images. Default: 40
-point_pixel
turn off solid-angle correction for square flat pixels
-printout
print pixel values out to the screen
-noprogress
turn off the progress meter
-nointerpolate
turn off the tricubic interpolation
-interpolate
turn on the tricubic interpolation, even for crystals
-round_xtal
use ellipsoidal crystal shape for spot shape calculation (approximate)
-square_xtal
use paralelpiped crystal shape for spot shape calculation (exact)
-binary_spots
cut off every spot at the FWHM, even intensity inside. not the default
-seed
manually set the random number seed. Default:
-mosaic_seed
different random number seed for mosaic domain generation. Default: