# Coordinate transformation¶

## Notation¶

$$X$$, $$Y$$ are square matrices with basis vectors as rows

$$X$$ … old, shape: (3,3) .. or (M,M) in general
$$Y$$ … new, shape: (3,3)
$$I$$ … identity matrix, basis vecs of cartesian system, shape: (3,3)
$$C$$ … transformation matrix, shape(3,3)
$$v_X$$ … row vector v in basis X, shape: (1,3)
$$v_Y$$ … row vector v in basis Y, shape: (1,3)
$$v_I$$ … row vector v in basis I, shape: (1,3)

## Row vs. column form¶

In the math literature you find a column oriented form, where matrix $$A$$ has the basis vectors as columns and $$x$$, $$b$$ are column vectors (M,1), such that a linear system is written as $$A x = b$$.

We use a row oriented form of all relations such that $$x A = b$$, where $$x$$ are $$b$$ are row vectors (1,M) and $$A$$ has basis vectors as rows. This is optimal for direct translation to numpy code, where we can use broadcasting.

All formulas here translate by $$(A x)^\top = x^\top A^\top$$ and $$(A^{-1})^\top = (A^\top)^{-1}$$.

Tools like linalg.solve and Lapack solvers assume the column oriented form, so you need to use linalg.solve(A.T, b.T).

## Theory¶

We have

\begin{align}\begin{aligned}v_Y Y = v_X X = v_I I = v_I\\v_Y = v_X X Y^{-1} = v_X C\end{aligned}\end{align}

where we call $$C = X Y^{-1}$$ the transformation matrix.

Every product $$v_X X; v_Y Y; v_I I$$ is actually an expansion of $$v_{X,Y,...}$$ in the basis vectors contained in $$X,Y,...$$ . If the dot product is computed, we always get $$v$$ in cartesian coords.

### Notes for the special case fractional <-> cartesian¶

With $$v_I$$ cartesian and $$v_Y$$ fractional coordinates, the transform fractional -> cartesian is the dot product

$v_Y Y = v_I$

from above. Cartesian -> fractional is

$v_Y = v_X Y^{-1}$

which is the solution of the linear system $$v_Y Y = v_I$$. It cannot get more simple.

## Implementation¶

We now switch to code examples, where $$v_X$$ == v_X.

In general, we don’t have one vector v_X but an array R_X of shape (N,M) of row vectors:

R_X = [[--- v_X0 ---],
[--- v_X1 ---],
...
[-- v_XN-1 --]]


We want to use fast numpy array broadcasting to transform all the v_X vectors at once. The shape of R_X doesn’t matter, as long as the last dimension matches the dimensions of C, for example R_X: (N,M,3), C: (3,3), dot(R_X,C): (N,M,3)).

Examples:

1d: R_X.shape = (3,):

>>> # R_X == v_X = [x,y,z]
>>> # R_Y == v_Y = [x',y',z']
>>> X=rand(3,3); R_X=rand(3); Y=rand(3,3); C=dot(X, inv(Y))
>>> R_Y1=dot(R_X, C)
>>> R_Y2=dot(dot(R_X, X), inv(Y))
>>> R_Y3=linalg.solve(Y.T, dot(R_X, X).T)
>>> assert allclose(R_Y1, R_Y2)
>>> assert allclose(R_Y1, R_Y3)


2d: R_X.shape = (N,3):

>>> # Array of coords of N atoms, R_X[i,:] = coord of i-th atom. The dot
>>> # product is broadcast along the first axis of R_X (i.e. *each* row of R_X is
>>> # dot()'ed with C)::
>>> #
>>> # R_X = [[x0,       y0,     z0],
>>> #        [x1,       y1,     z1],
>>> #         ...
>>> #        [x(N-1),   y(N-1), z(N-1)]]
>>> # R_Y = [[x0',     y0',     z0'],
>>> #        [x1',     y1',     z1'],
>>> #         ...
>>> #        [x(N-1)', y(N-1)', z(N-1)']]
>>> #
>>> X=rand(3,3); R_X=rand(5,3); Y=rand(3,3); C=dot(X, inv(Y))
>>> R_Y1=dot(R_X, C)
>>> R_Y2=dot(dot(R_X, X), inv(Y))
>>> R_Y3=linalg.solve(Y.T, dot(R_X, X).T).T
>>> assert allclose(R_Y1, R_Y2)
>>> assert allclose(R_Y1, R_Y3)


Here we used the fact that linalg.solve can solve for many rhs’s at the same time (Ax=b, A:(M,M), b:(M,N) where the rhs’s are the columns of b). The result from linalg.solve has the same shape as b: (M,N), i.e. each result vector is a column. That’s why we need the last transpose.

3d: R_X.shape = (nstep, natoms, 3):

>>> # R_X[istep, iatom,:] is the shape (3,) vec of coords for atom iatom at
>>> # time step istep.
>>> # R_X[istep,...] is a (nstep,3) array for this time step. Then we can use
>>> # the methods for the 2d array above.
>>> X=rand(3,3); R_X=rand(100,5,3); Y=rand(3,3); C=dot(X, inv(Y))
>>> R_Y1=empty_like(R_X)
>>> R_Y2=empty_like(R_X)
>>> R_Y3=empty_like(R_X)
>>> for istep in range(R_X.shape):
>>>     R_Y1[istep,...] = dot(R_X[istep,...], C)
>>>     R_Y2[istep,...] = dot(dot(R_X[istep,...], X), inv(Y))
>>>     R_Y3[istep,...] = linalg.solve(Y.T, dot(R_X[istep,...], X).T).T
>>> assert allclose(R_Y1, R_Y2)
>>> assert allclose(R_Y1, R_Y3)


Here, linalg.solve cannot be used b/c R_X is a 3d array and not a matrix. Direct dot products (R_Y1 and R_Y2) involve calculating the inverse, which is unpleasant. Hence, the numerically most stable method is to loop over nstep and use the 2d array method using linalg.solve (R_Y3).

The above loops are implemented in flib.f90 for the special case fractional <-> cartesian, using only dot products and linear system solvers from Lapack in each loop.