I'm a chemistry student and I'm currently trying to write a Python script that helps to solve a crystal structure based on experimental values.
The equation for solving a cubic crystal (the simplest case) is:
$$ 4 \cdot \frac{ \sin^{2}{\theta} }{\lambda^2} = \frac{1}{d^2}=\frac{1}{a^2}\cdot\left(h^2 + k^2 + l^2\right) = P\cdot(N)$$
The experimental result of an x-ray-diffraction gives $2\theta$. Given this, it's pretty simple to calculate $1/d^2$, as $\lambda$ is known too. That's the wavelength of the x-ray beam.
The unknown values of chemical interest are $a, h, k$ and $l$. We know that there is a constant $P$ which leads to the $1/d^2$ values by multiplication with a natural number $N$. Furthermore, it must be possible to represent that natural number by a sum of squares: $\left(h^2 + k^2 + l^2\right)$.
Given this, I thought that I have to determine the greatest common divisor for all the $1/d^2$ values. Unfortunately, these numbers are not natural numbers. See Table 1 for an example of values.
My problem is, that I don't have a good idea to do so. My attempt was to divide each $1/d^2$ value by a natural number $x$ and save the result to a list. Then divide by $x+1$ and so on. The limit would be set arbitrarily. After dividing the first $1/d^2$ value by that range of numbers, I'll move on to the second one, and so on.
After having done this for all the $1/d^2$ values I'll check which value is the most common one and suppose that that is my greatest common divisor.
Yes, I am totally aware of the fact, that you are crying now because that method is very bad. There are several problems with this and I know that, but it was the only thing I came up with. And that's why I'm here, I hope that anyone can help me with this question:
How do I compute the greatest common divisor for my $1/d^2$ values in a mathematically correct manner?
Thanks in andvance!
Update
I thought about the problem once more and I want to explain in greater detail what the Python script should be able to do. With these information, it's also more simple to help I guess.
Transform the $2\theta$ values into $1/d^2$ values. (solved)
Search the greatest common divisor for these values.
Check if $N$ results in natural numbers that can be produced by a sum of three squared natural numbers $\left(h^2 + k^2 + l^2\right)$. A mandatory condition for the solution of the crystal structure.
If this is the case, save the values for $h, k$ and $l$ for each $1/d^2$ value for future printing.
If this is not the case, take the second greatest common divisor and test again.
Calculate the value for $a$. (solved)
Taken out the fact, that I'm still not able to calculate the greatest common divisor, I think that my plan could work out well. I wanted to set up a simple algorithm that populates a list by performing the $\left(h^2 + k^2 + l^2\right)$ calculation and storing the used numbers as well as the result. Then, after finding a GCD and computing $N$ for each value, I'll check if there is any value in the list of $N$ that is not represented in the list of the $\left(h^2 + k^2 + l^2\right)$ results. With this, I know if I have to take the second greatest common divisor or not.
Again, I'm happy for every shared idea on how to get the greatest, or/and even the second or third greatest common divisor. And please keep in mind, that I'm not a mathematician. ;-)
Table 1
\begin{array} {|r|r|} \hline No & 2\theta & 1/d^2 \\ \hline 1 & 15.845 & 0.24232239557963048 \\ 2 & 22.480 & 0.4846105998701098 \\ 3 & 27.623 & 0.7269290249651083 \\ 4 & 32.003 & 0.9692579815933621 \\ 5 & 35.901 & 1.2115343992719871 \\ 6 & 39.463 & 1.4538572457953847 \\ \hline \end{array}
For the given data, you can determine the smallest common divisor simply by inspection: all six entries are extremely close to the arithmetic progression of multiples of the first entry.
For other data, you can perform the Euclidean algorithm (or a multi-argument variant), but the main point is that if you ever reach a sudden, dramatic drop in the size of a number, then it's likely that you have computed what should be zero. You can then spin off a new problem where you replace the extremely small number with zero, and see if you get a result that looks good.
(if you don't, you might go back to the original problem and continue without doing the substitution)
Another mildly ad-hoc approach is to compute the ratio of every term to the first, and use continued fractions to estimate the proportion with a rational number. Hopefully these are the true ratios, at which point you can use an exact GCD algorithm.
The systematic (and probably most reliable) approach to the problem is by lattice reduction; e.g. the basic approach is to set up a matrix
$$\begin{matrix} W & 0 & 0 & 0 & 0 & 0 & 0.24232239557963048 \\ 0 & W & 0 & 0 & 0 & 0 & 0.4846105998701098 \\ 0 & 0 & W & 0 & 0 & 0 & 0.7269290249651083 \\ 0 & 0 & 0 & W & 0 & 0 & 0.9692579815933621 \\ 0 & 0 & 0 & 0 & W & 0 & 1.2115343992719871 \\ 0 & 0 & 0 & 0 & 0 & W & 1.4538572457953847 \\ \end{matrix}$$
and use a lattice reduction algorithm to reduce the rows of the matrix. Hopefully, five of the rows will have an extremely small value in the last column, and the sixth row will have in its last column an estimate of the greatest common divisor of your values.
Setting the value of $W$ is an art. You want:
The second bullet point is the reason for the extra columns, by the way — they keep track of how much arithmetic you've done on the rows, and thus estimate how much error accumulation has happened.
People who actually know the subject can do math to inform their value of $W$. Us mortals should just run lattice reduction a lot and keep decreasing $W$ until they get a result that looks right, and stop when the the results of lattice reduction start becoming garbage.