Interfacing ase and 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!