# 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.