Coordinate transformation¶
Notation¶
\(X\), \(Y\) are square matrices with basis vectors as rows
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
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
from above. Cartesian -> fractional is
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.