# OlexSys : Twinning

--------------------------

Comment by Colin Please:

Given an observation $F_{obs}(hkl)$, we are seeking to find a calculated distribution that fits this data. If there is no twinning then we can seek
$$|F_{obs}(hkl)|^2 \approx |F_{calc}(hkl)|^2.$$
However, if there is a single region of twinning then we might seek to fit the data using
$$|F_{obs}(hkl)|^2 \approx |F_{calc}(hkl)|^2 + \beta |F_{calc}(R*hkl)|^2$$
where the scalar $\beta$ and the rotation matrix $R$ are additional fitting parameters.

A possible route forward it to use the existing procedures to do the fitting assuming that no twinning has occured (twinning is typically a perturbation to the system). This gives us a possible $F_{calc}^n(hkl)$.

The major difficulty with finding $R$ is that the residual is a very sensitive function of the angle. Hence, a gradient method to find the optimum will proabably not work. An alterntive is to explore a small set of $R$ that are expected to be likely. The current method appears to do this by looking at flipping and other simple rotations.

An alternative approach might be to use the estimated $|F_{calc}^n(hkl)|$ and subract this from the complete original data set $|F_{obs}(hkl)|$ to find a residual. We are now seeking to find a rotation matrix $R$ and a scaling factor $\beta$ that fits this data. We might note that the fitted data $|F_{calc}^n(hkl)|$ consists of a number of intensities at the integrer points of $(hkl)$. Hence we coudl look at the residual and order the intensity of the points in this data (note that many of these may not be at integer points). We can then compare this list with the ordered list in the calculated values and see if the distance of the point from the origin is the same. From this we can determine a small list of possible candidate points in the residual that might correspond to a point in the calculated data. Comparing the position of these candidate points with the calculated points then gives a short list of candidate $R$ to explore (as well as rough estimates of the possible $\beta$). Given this short list local refinement can occur t see if the fit is good.

Note the candidate list could readily identify several possible $R$ and it might be necessary to consider more than one region of twinnning but the procedure above would seem to allow this.

--------------------------

## Tuesday morning discussion

A fundamental issue is that we only have the hkl file; it isn't feasible to have all of the intensities (including all the points in space), we will only have the intensities at integer hkl file. As a result, we can't use the full data

Therefore the interesting issue is what rotation matrices will rotate some proportion of the lattice points onto other lattice points.

We are looking for possible overlaps robably within 5 lattice points, definitely within 10 lattice points..

## Tuesday afternoon thoughts

* Annuli -- given a tolerance $\epsilon$ and a given basis for a crystal, what are the radii which give nontrivial number of rotations on to the top of itself.

* Amit considered a triplet of peaks in $R^3$ with the property that the first peak is mapped onto the second and the second onto the third, and he then thinks that this allows to find the axis and the rotation.

* Michael: The other idea was to take the residual $Residual = I_{obs} - I_{calc}$ and then to consider this as a subset of the rotated lattice and try to find the full rotated lattice from this info or, at least, certain rotations who will allow this, that is find $R_\theta$ such that $Residual \subset R_\theta \cdot Gamma$.

* Could build on the annuli work since this can suggest possible ways of making twins in synthetic data

* A more theoretical problem: Given a 2--dimensional lattice $\Gamma=\mathbb{Z}a+\mathbb{Z}b$, and a rotation $R_\theta\in SO(2)$, under what condititions is $R_\theta\Gamma\cap \Gamma$ nontrivial? How big can the intersection be? Note that for a given $x,y\in \Gamma$, $y=R_{\theta} x$ for some $\theta$ if and only if $|x|=|y|$, which helps restrict the search. In practise should work with a tolerance.

* Look into Diophantine approximation

* Do a proper literature review and discover that this has all been known since the 1930s.

* Need to read the Rotax paper

https://onlinelibrary.wiley.com/doi/full/10.1107/S0021889802000249?sentby=iucr

* Paper relating Statistics and Twins:

https://journals.iucr.org/d/issues/2006/01/00/ba5089/index.html

* Maybe another useful paper on twinning:

https://core.ac.uk/download/pdf/205704459.pdf

## Wednesday Thoughts

* Cameron
* 2D: For a given $b$, what is the minimum radius for which non-trivial twinning (with tolerance) occurs? Blue figure
* 2D: From this, which $b$ lead to non-trivial twinning (with tolerance) within say 5 $h$ coordinates. Blue and yellow figure.
* James
* 2D: The $b$ vectors which can lead to non-trivial twinning (tolerance=0) lie on circlines with formulae
$$
b_{1}=\frac{-\alpha}{2\beta},
$$
$$
\gamma\left(\left(b_{1}+\frac{\beta}{\gamma}\right)^{2}+b_{2}^{2}\right)+\alpha-\frac{\beta^{2}}{\gamma}=0,
$$
where $\alpha,\beta,\gamma\in\mathbb{Z}$. This comes from taking two lattice points and setting their distances from the origin to be equal.
*
* Michael:
* In 2D, Cameron's code creates a bunch of fake data for the XRD data, giving a set of points in Euclidean space of intensities, with or without twinning. We run this including twinning to get an approximation for $F_\text{obs}$.
* Running this code without twinning gives us a result for $F_\text{calc}$, similar to what Olex2 would predict if it was ignoring twinning. We identify possible twinning cases by creating a set of residuals equal to $F_\text{obs} - F_\text{calc}$.
* These will lie on lattice points of the original lattice, and the strength of the intensity is important for determining the value of $\beta$ later.
* Our first aim is to determine the twinning angle from this. We start with the original lattice (the $F_\text{calc}$), and then start to rotate it about the centre point. For each angle, we calculate how many of the residual points line up with this rotated lattice (up to a tolerance).
* We continue to rotate up to a maximum of 180 degrees. Once we do, we can plot the number of aligned points as a function of the angle
* We see that the number of aligned points is almost everywhere zero, with occasional spikes where we get sudden alignment. Even sampling 1000 points between 0 and 180 degrees, this still happens and shows how non-smooth this problem is
* Also of note is the symmetries of the lattice also cause overlapping, to exclude this we need to know the symmetries a-priori in order to exclude these [name=jph this narrows down your search window?]
* It should then be simple to identify which of the peaks are truly due to twinning, hence identifying their angle
* We should be able to identify the value of $\beta$ from these residual points. Recall that
$|F_\text{obs}|^2 = |F_\text{calc}|^2 + \beta |R F_\text{calc}|^2$,
where $R$ is the rotation matrix. Knowing the value of R means we can determine the value of $R F_\text{calc}$ at each of the residual points, which gives us a set of points where we can solve the equation above to find $\beta$. Doing this for all of the points, we can find out what the value of $\beta$ is with an identifiable error.

* Amit: In a group with Luke, James and Christian we have discussed the following "algorithm" (we assume that we know the original lattice structure):
**Step 1:** Create a list that assigns the distance to the origiin of any grid point to the point. We will commonly call these distances "Radii"
**Step 2:** For a given possible radius, choose a tolerance parameter $\epsilon$ and consider all the grid points that appear in the annulus with internal radius $R-\epsilon$ and external radius $R+\epsilon$ (we can choose different upper and lower tolarence parameters). Denote this set by $G_{R,\epsilon}$.
**Step 3:** If $\left| G_{R,\epsilon}\right|>1$ then there exists two a lattice point which can be rotated to another lattice point. Here we have two options:
*Option 1:* We are in $2D$. We create a list of pairs of points in $G_{R,\epsilon}$ together with the cosine of their angles (computed using the standard inner product). All pairs that have the same cosine of the angle between them (as long as we are restricted to half a lattice, I believe) are those who will contribute to twinning. Moreover (if I am not mistaken - Amit), the number of such pairs is exactly the number of intensified points we'll discover.
*Option 2:* We in $3D$. In this case we can continue in a similar way to the above: We create a list of pairs of points in $G_{R,\epsilon}$ together with the cosine of their angles. This time, however, we can use this angle to know exactly how many points are mapped to others with the same rotation. However, any pair in $G_{R,\epsilon}$ yields a possible twining. Creating the plane that passes through these points and the origin, we use the angle between these vector to craete a rotation around the normal of this plane with the same angle.

* About angular tolerance: In the above we allowed the pairs to not necessarily be on the same radius (radial tolerance). Howver, we can also allow the grid points to not be perfectly aligned. This corresponds to choosing rotation matriced with a $\pm \delta$ error to the exact angle computed in the previous procedure.

## Thursday Thoughts

* How many rotations are there between two points on a given sphere? A rotation $R\in SO(3)$ can be parametrised by an axis $n$ (a unit vector) and a rotation angle $\theta$. The explicit rotation matrix can be given as
\begin{equation}
R_{ij}(n,\theta)=n_in_j+(\delta_{ij}-n_in_j)\cos\theta-\epsilon_{ijk}n_k\sin\theta.
\end{equation}This parametrisation is one way to realise the isomorphism $SO(3)\simeq \mathbb{RP}^3$. Given $u,v$ on the unit sphere $S^2$, there is a one parameter family of rotations joining them. The axis is
\begin{equation}
n(t)=\frac{u+v}{2\cos(\alpha/2)}\cos t +\frac{u\times v}{\sin\alpha}\sin t, \hspace{1cm} t\in [0,\pi],
\end{equation}where $\alpha=\cos^{-1}(u\cdot v)$ is the angle between $u$ and $v$. The angle was found by Cameron to be
\begin{equation}
\cos\theta(t)=\frac{\sin^2 t\cos^2(\alpha/2)-\sin^2(\alpha/2)}{\sin^2 t\cos^2(\alpha/2)+\sin^2(\alpha/2)}.
\end{equation}The minimum value $\theta=\alpha$ occurs at $t=\pi/2$ and is exactly the angle $\alpha$. The maximum value is $\theta=\pi$ and occurs at $t=0,\pi$. $(n(t),\theta(t))$ somehow parametrises an ellipse when viewed in terms of quaternions (see below). This suggests the following algorithm for identifying twin laws for a real crystal (which as I understand has basically already been tried by Norbert and Laura):

* Calculate residual lattice $F_{\text{obs}}-F_{\text{calc}}$, identify largest vector $a$. Calculate length $r$ of $a$.
* Identify points $\{b_1,b_2,\dots,b_N\}$ in the lattice with the same radius $r$ (up to a tolerance).
* Calculate the 'rotation ellipses' $E_1(t),E_2(t),\dots, E_N(t)$, where $E_i(t)$ represents the $S^1$ of rotations mapping $a$ to $b_i$. In total there are $N$ circles of rotations.
* Test each of these rotations as a possible twin law, for example by investigating the fate of the top 30 elements of the residual lattice.

There are a few problems with this but it's a decent start...

* Sets of Rotations can be written as ellipses representing quaternians, where the direction is the axis and the angle is given by the distance (somehow). These ellipses have major axes touching the unit sphere and minor axes within the unit sphere. If such ellipses were 'close' at a point, that point would represent a possible rotation for both sets. If we were to gain 'all' ellipses for our strongest reflections and look for points within the unit sphere which have 'many' close by ellipses, those points would represent very likely twin laws!