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[0]):
>>>     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.