How to cut a slice in a 3D grid

Let define a 3D-grid, the pixels are not necessarily cubic. This grid contains the electron density resulting from a fourier transform. On pixel define a basis A(a1,a2,a3) of the grid. the length of the grid cover the unit cell in 3 directions.

The task is to cut a 2D-grid plane in this 3D object. Then, a contour plot can be drawn using the plane. A figure showing the transform between two 2D grids shows the dificulty of the task.

We need first to define the plane with 3 points P, Q and R (atoms) lying in this plane. The first thing to do is to choose an origin and define an orthonormal basis S(s1,s2,s3) for it. Let be P the origin of the basis. The A basis is not orthonormal, the ortohonormalisation matrix O is first use to get the coordinates of P, Q and R in an orthonormal coordinate system AO (ao1,ao2,ao3). The reason is the use of the cross product to define the basis of the plane.

The basis (s1,s2,s3) of the plane is define as:
\(\overrightarrow{s2}=(\overrightarrow{Q}-\overrightarrow{P})\times (\overrightarrow{R}-\overrightarrow{P}) \)
\(\overrightarrow{s3}=\overrightarrow{s1} \times \overrightarrow{s2} \)

Now, we can define a tensor than going to transform orthonormal crystal coordinates into the new basis of the plane. This tensor operates on from the basis AO to S.

\( T = \begin{bmatrix}
s1_x^{ao} & s2_x^{ao} & s3_x^{ao} \\
s1_y^{ao} & s2_y^{ao} & s3_y^{ao} \\
s1_z^{ao} & s2_z^{ao} & s3_z^{ao} \\
\end{bmatrix} \)

The two tensors T and O are going to map the pixels from the 3D-grid to the 2D-grid. The next step is to define the grid fot the 2D-plane. Let define a plane width a width w of 16 ang and a ratio width/height r of 16/10. The resolution of the grid will be 0.1 ang, meaning a 160×100 grid plane. The pixel (i,j) will have the coordinates (x,y,0) in S with:

\( x = i \cdot w/(160/\|\overrightarrow{s1^{ao}}\| ) \)
\( y = j \cdot (w/r)/(100/\|\overrightarrow{s2^{ao}}\| ) \)

Now we map the coordinates of each (i,j) pixel plane to the coordinates (u,v,w) in the 3D-grid.

\( (u,v,w) = \mathbf{O}^{-1} ( ( \mathbf{T}^{-t} \begin{bmatrix}x \\ y \\ 0 \end{bmatrix} ) + \overrightarrow{P^{ao}} ) \)

This was the easy part. The (u,v,w) coordinates need to be matched to a pixel in the 3D-grid. To go from fractional coordinates to the pixel index, it is just necessary to multiply by the length in grid size in each direction. The problem is that the result is not an integer and an interpolation is needed. The easiest interpolation is to the nearest pixel, practically, we round the result (u,v,w) to the nearest integer. An other easy solution is to use tri-linear interpolation. Platon is using a more complex quadratic or cubic interpolation.

The resolution of the grid plane is arbitrary and have to be chosen carefully according to the original grid resolution. In case of wide angle and a cut through the diagonal, there is a loss of resolution. For example with a 0.3 ang grid resolution in the 3D-grid, the maximum resolution will be about 0.7 ang with a diagonal cut in the trigonal space group. 0.3 ang is the resolution used by platon. Hence the use of heavy interpolation in order to smooth the contour plot. Using a linear interpolation, a minimum resolution of 0.1 ang is necessary, 0.05 is perfect. It represents a 3D-grid of about 1M-8M pixels. The fftw library for the fourier transform handle this quite easily (0.2s-1.5s).

Interpolation to the nearest integer of a 0.3 ang 3D-grid :

High resolution 3D-grid of 0.05 ang. The 2D-grid resolution is kept at 0.01 ang to avoid losing any information and a linear interpolation has been used.

Comments are closed.