HKL reflections for X-Ray Crystallography

[Phonon]

I have a feeling that there aren't many X-ray crystallographers here. I wrote this for a 4th year course in x-ray crystallography. At the moment, it only works for orthorhombic crystals. You input the lengths of the cell axes in Angstroms and the program will output the distance between reflection planes (hkl) and the reflection angle theta, for h, k, l between and including -10 to 10.

I'm in the middle of writing another one that will output a Patterson function for cracking the phase problem in x-ray crystallography.

code:

cls function arcsind (r : real) : real     result arctand (r / sqrt (1 - r ** 2)) end arcsind var wavelength : real var a : real var b : real var c : real var theta : real var nMax : int := 0 var sinT : real var nUse : int := 1 var n : int := 1 var stremout : int open: stremout , "data.txt" , put put "Enter the wavelenth of radiation in Angstroms: " get wavelength put "Enter the length of the unit cell axis a in Angstroms: " get a put "Enter the length of the unit cell axis b in Angstroms: " get b put "Enter the wavelength of the unit cell axis c in Angstroms: " get c put " " put "The various reflections will written to data.txt." put " " put "The output format is: h, k, l, d(hkl), theta (in degrees), order of reflection." for h : -10 .. 10     for k : -10 .. 10         for l : -10 .. -1             var dist : real             dist := (1/sqrt(((h*h)/(a*a))+((k*k)/(b*b))+((l*l)/(c*c))))             n := 1             loop                 sinT := ((n*wavelength)/(2*dist))                 exit when sinT > 1                 theta := arcsind (sinT)                 put: stremout, h : 3 , k : 3 , l : 3 , dist : 15 , theta : 15 , n : 7                 n := n + 1             end loop         end for                    end for                 end for   for h : -10 .. 10     for k : -10 .. 10         for l : 1 .. 10             var dist : real             dist := (1/sqrt(((h*h)/(a*a))+((k*k)/(b*b))+((l*l)/(c*c))))             n := 1             loop                 sinT := ((n*wavelength)/(2*dist))                 exit when sinT > 1                 theta := arcsind (sinT)                 put: stremout, h : 3 , k : 3 , l : 3 , dist : 15 , theta : 15, n : 7                 n := n + 1             end loop         end for                    end for  end for close : stremout

[codemage]

You've lost me. I don't have the free time I used to for crystallography.

[Delos]

Next step - visualization! And then from there, why not interpret the pattern too? Who knows where that will take you. Proteins at a glance? I'm sure that would be terrible fun.
Interesting work to say the least. As you said, not many of us are Crystalographers, and fewer still would be able to understand the formulae you used or more specifically, be able to appreciate its purposes. For those interested, why don't you give some sample output and explain what the little numbers mean, just so we know what we're dealing with when we see "unit cell axis a".

[Phonon]

The Patterson function program which I write will give graphical output and I'll post it here. Here's an example of what this program does though:

For all of you non-crystallographers, the unit cell is the simplest unit of repeat that a crystal can have. Crystals are made up of repeating units of atoms or molecules or ions (or some combination thereof - and of course, molecules are made up of atoms, but I'm referring to nonbonded atoms). In an orthorhombic cell, which is what this program works for, all angles of the unit cell are 90 degress, but the lengths of each cell axis can be different, so that's what you input. The radiation wavelength is also input and generally they're around 1 angstrom. An angrstrom is 10^-10 m, or 0.1 nm.

So, here's sample output data for a = 1, b = 2, c = 1.5, wavelength = .71073. Not all of it is shown.

============================
2 0 1 0.474342 48.518982 1
2 0 2 0.416025 58.670491 1
2 1 1 0.461538 50.350073 1
2 1 2 0.407307 60.747602 1
2 2 1 0.428571 56.015109 1
2 2 2 0.384111 67.692966 1
2 3 1 0.386494 66.846982 1
=============================

The first three numbers are h,k and l. Those numbers make up what are known as the Miller indices. The Miller index is a way of denoting a given plane of reflection. If you imagine a unit cell which is a cube, with an atom at each of its 8 corners, then you can see that there you can draw a plane through the atoms at the four corners of x=0 (the atomic coordinates are (0, 0, 0), (0, 1, 0), (0, 1, 1), and (0, 0, 1).) There is an identical plane at x=1. The Miller index for this plane has h = 1, k = 0, and l = 0, or (hkl) = (100). This is because the reciprocal of the distance between which the plane intersects the x-axis is 1 (it intersects 1 per axis). Since the planes are parallel with the y and z axes, then these planes take an infinity until they intersect these axes again. The reciprocal of infinity is 0, so for the y-axis k = 0 and for the z-axis l=0.

The next number in the output is the distance between each (hkl) plane. As you would expect, numbers with larger h, k and l values have smaller distances between them. The next value is the kind of thing you'd need to search the net for to get a visual of ("Bragg's law"), but it is the reflection angle in degrees. For higher values of h, k and l, the reflection angle must get correspondingly smaller and if you find an image showing Bragg's law, you'll understand why. I haven't included this in the program, but higher angle reflections give weaker intensities, so I could even include a prediction of expected intensity and generate a map of the reflections in reciprocal space like one of the HKL tools in the WinGX program (which is what Delos may be suggesting). The last number is the order of the reflection and it has to do with how many times the wavelength can fit in between the (hkl) planes before reflecting. Again, if you find an image showing Bragg's law, you'll see why higher angles allow for higher values of n (the order of reflection).

Delos, what kind of visualization were you thinking of? A reciprocal space map showing predicted intensities of each reflection in reciprocal space, like in WinGX? For that, I'd probably have to make a substantial modification where I actually inputted atom types, atomic positions, calculated structure factors, and so forth. It's a cool idea though. I could probably get the program to generate .HKL files which I could then plug into WinGX and solve!

[Delos]

[nervous laughter]

The closest I've ever come to thinking about crystallography theoretically is some basic NMR (H and C) and IR theory. Yeah, that far away...hence as my knowledge is so utterly primitive, I was simply referring to one of those nifty diagrams that look like the Polka-dot Monster exploded on them. (c.f. the oh-so famous Franklind-Watson-Crick structure of DNA). The type that people sit in front of for hours and then in a fit of monumental insight, divine a structure from.

As you can see, I'm way off. My brother is actually doing a lot of this sort of stuff at UofT right now, and is also taking a Crystallography course. He'd appreciate this stuff way more!

[Phonon]

You mean something like this?

[Delos]

Looks purdy. Now what are we looking at there?

[Phonon]

Reciprocal space reflections.

So this is one of the (hkl) planes, where you're viewing along the k = 0 plane. The different h values are off the the left and right directions and the various l values are up and down. So, for exmple, the dot on the first column to the right and the second row down is (hkl) = (102). While the dot on the third column to the right and the third column down is (hkl) =(303). So this is a map of the intensities from the reflections of the various HKL planes. My program tells you the distances between HKL planes and the reflection angle, which isn't shown here. This picture doesn't show distances or angles, but it shows the intensities of the reflections and the legend can be seen on the left side, with larger numbers (and dots) meaning higher intensity. Since the distance between planes decreases with increasing values of h, k and/or l, the reflection angle increases (look at Bragg's law). As I pointed out, the reflection intensity decreases with increasing angle, thus we would expect a reflection like (10 0 8) to be of low intensity.

So, take a look at the 10th column to the right and the 8th column down on the map. I have the program window open and I can click on the corresponding dot to get the intensity, but if you look at the legend you can see that the reflection is of very low intensity. It's a blue "+", so the relative intensity is less than 10, while the most intense reflections here are between 2500 and 4999 (when I click on the intense ones they're about 3000). The "reflection" at (10 0 8) is probably just experimental error and the crystal isn't reflecting from that plane.

So, what does all of this tell you? Different crystal systems will have symmetries which have reflections which cancel each other out due to the electromagnetic waves having opposite phases (destructive interference). These absences are known as "systematic absences" and they give you clues as to the symmetry of the crystal. The WinGX program looks at the reflection data and looks for various systematic absences to see whether or not it can deduce certain symmetry elements, such as a glide plane or screw axis. If the program can solve the total symmetry of your crystal then it can use the reflection data to produce a map of the electron density in your crystal (after solving what's known as "the phase problem" which I won't get into now). This is where the chemist's skills come in. The chemist looks at the electron denisty map and starts assigning a structure to it. You can then use the program to refine the structure you've assigned against the electron density map. Basically, the process is iterative. You assign atoms to a spots of electron density based on what you know about chemisty (for example, you can easily pick out benzene rings) and have the program refine the electron density map based on your changes. You then get a better electron density map and keep updating it. You can calculate what a perfect set of reflections would look like based on your assigned structure and compare it with the observed reflections. The relative error between the two gives you an idea of how close your structure is to likely being correct. Generally, once you get the error to less than 5% your structure is satisfactorily solved.