2021-03-28

# Interfacing ase and i-pi

### Context

The Atomic Simulation Environment (`ase`) is a python package providing common classes and various tools for describing and simulating atomistic systems. `i-PI` is a python package in the same domain, but focused more on path-integral molecular dynamics. For the `gknet` project, I'd like to run some simulations with `i-pi`, but my backend is based on `ase`. This note contains some pointers towards interfacing those two. (I'd recommend also checking out the documentation of `i-pi`!)

## Coordinate system

For periodic systems, `i-pi` internally works in a "canonical" coordinate frame where the first lattice vector points along the x axis, the second one lies in the x-y plane, and the third vector is arbitrary. (This makes the unit cell matrix triangular, which is convenient.)

Orienting the system in such a way can be achieved in `ase` with:

``````def orient_atoms(atoms):
new = atoms.copy()
new.set_cell(atoms.cell.cellpar())
new.set_scaled_positions(atoms.get_scaled_positions())

return new
``````

This exploits the fact that `ase` uses the same convention as `i-pi` to reconstruct a unit cell from `cellpar` (the lengths and angles of the basis vectors).

## Notations for unit cell matrix

In `i-pi`, the basis vectors are the columns of a `cell` matrix, whereas in `ase`, the `cell` contains those vectors as rows. To convert from `ase`, one therefore has to use the transpose.

## Input format

`i-pi` supports the `.xyz` format for input geometries. However, `xyz` doesn't have a standard way to write down the basis vectors for periodic systems. In `i-pi`, there are three options to specify the basis: `H`, `abcABC`, and `GENH`.

For `H` and `abcABC`, `i-pi` expects the system to be rotated "canonically", for `GENH` this is waived, as the rotation is performed internally on import.

The `H` mode expects the unit cell as the `.flatten()`-ed version of the `i-pi`-style `cell` matrix, i.e.:

``````def get_cellh(atoms):
"""Get H-mode cell spec from a canonically-oriented atoms object."""
return atoms.get_cell().T.flatten()
``````

Note the tranpose to convert from `ase` row style. This format is also used in `input.xml` to specify the unit cell for barostats, so it's useful to use.

The `abcABC` mode expects the cell parameters:

``````def get_cellabc(atoms):
"""Get abcABC-mode cell spec from a canonically-oriented atoms object."""
return atoms.get_cell().cellpar()
``````

Since the cell parameters don't depend on the orientation of the structure, no alignment needs to be performed here.

Finally, the `GENH` mode expects the unit cell in row style (!). The `Atoms` object doesn't have to be re-oriented:

``````def get_cellgenh(atoms):
"""Get GENH-mode cell spec from an atoms object."""
return atoms.get_cell().flatten()

``````

However, since the re-orientation will be done anyway by `i-pi`, I'd recommend starting with a properly oriented system, and not using this mode. Anecdotally, it seems that the `ase` way results in less numerical noise in components of the positions that should be zero.

In all cases, the basis is represented as a comment line, formatted as:

``````# CELL(\$MODE): \$SOME \$NUMBERS cell{\$UNIT} positions{\$UNIT}
``````

In other words, to get from `ase.Atoms` to `i-pi`-style `xyz`, one needs to run one of:

``````
def xyz_cellh(filename, atoms):
rotated = orient_atoms(atoms)
comment = f"# CELL(H):     " + "     ".join([f"{x:.5f}" for x in get_cellh(rotated)])
comment += r" cell{angstrom}"
comment += r" positions{angstrom}"

write(filename, rotated, format="xyz", comment=comment)

def xyz_cellabc(filename, atoms):
rotated = orient_atoms(atoms)
comment = f"# CELL(abcABC):     " + "     ".join(
[f"{x:.5f}" for x in get_cellabc(rotated)]
)
comment += r" cell{angstrom}"
comment += r" positions{angstrom}"

write(filename, rotated, format="xyz", comment=comment)

def xyz_cellgenh(filename, atoms):
comment = f"# CELL(GENH):     " + "     ".join(
[f"{x:.5f}" for x in get_cellgenh(atoms)]
)
comment += r" cell{angstrom}"
comment += r" positions{angstrom}"

write(filename, atoms, format="xyz", comment=comment)

``````

## Small gotchas

### Logging

The positions are written down canonically oriented in log files and are also wrapped back into the unit cell. As far as I can tell, the positions sent to calculators are not wrapped, so if one wants to compare for some reason, running `atoms.wrap()` makes positions (mostly) comparable. There might still be some differences close to the unit cell boundaries, depending on tolerance settings.

### Stress with `SocketClient`

If the stress is required, for instance for NPT, the client needs to be started with `client.run(atoms, use_stress=True)`.

Please note: this document is very much a work-in-progress, and should be taken as such. If you find anything wrong just drop me a line at `mail@marcel.science`, or let me know on twitter. Same if you find anything wonky with the website, it's all new!