3
$\begingroup$

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.

  1. Transform the $2\theta$ values into $1/d^2$ values. (solved)

  2. Search the greatest common divisor for these values.

  3. 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.

  4. 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}

$\endgroup$
4
  • $\begingroup$ Your experiment cannot produce 17 significant digits. How many can you trust? Moreover, unless you can say something about P, the problem is ill-defined: you can choose arbitrary N (expressible as a sum of three squares) and calculate a corresponding P. $\endgroup$ Commented Nov 8, 2017 at 0:03
  • $\begingroup$ You are absolutely correct. I was just pasting the calculated results in here, without taking care of the significance. But as you are asking, I can only trust 5 digits, which are the five significant digits of the $2\theta$ values. Well, $P$ equals $1/a^2$ and thus is directly correlated to $a$, which is the base length of an elementary cell. So I cannot say much about $a$, besides that it should be a real number in de range of about $5 Å$. $\endgroup$ Commented Nov 8, 2017 at 0:10
  • 1
    $\begingroup$ @Nick: Besides this, isn't it correct to calculate with the exact values as long as possible and only reduce the end result to the correct amount of significant digits? At least that's what we're told in chemistry. And the provided values are only a between result to show you these real numbers I'm dealing with. $\endgroup$ Commented Nov 8, 2017 at 0:21
  • $\begingroup$ P.S. you're looking for the greatest common divisor, not the smallest. The result is 'small' because its a common divisor, but you only want the result to be as small as necessary. Thus, you want the greatest among all common divisors. $\endgroup$ Commented Nov 8, 2017 at 0:51

1 Answer 1

Reset to default
1
$\begingroup$

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:

  • $W$ to be small enough so that lattice reduction will do enough work to find the answer.
  • $W$ to be large enough that lattice reduction will stop doing work before errors accumulate enough to spoil the result

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.

$\endgroup$
1
  • $\begingroup$ Thank you for your detailed answer. Unfortunately, as a chemist, most of you answer is like chinese for me. :-/ Could you please explain what lattice reduction means and how this relates to my problem? That would be really great. $\endgroup$ Commented Nov 8, 2017 at 9:41

You must log in to answer this question.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.