atoms out of simulation cell when converting to ase.atoms object

[naveenmohandas]

Hi, I was trying to convert a dislocation to an ase.atoms object. The conversion works fine however, in the ase.atoms object the atoms lie outside of the simulation box. Is there a proper way to move the atoms into the simulation box when converting dislocations created in atomman? and also when converting it back?

# code snippet after creating monopole dislocation as per https://www.ctcms.nist.gov/potentials/atomman/tutorial/4.9._Dislocation_configurations_generator.html

aseatoms, prop = disl_system.dump('ase_Atoms', return_prop=True)

from ase.visualize import view
view(aseatoms)

[lmhale99]

Hi,

I just uploaded a patch to this repository to fix this. The patch is also shown below if you want to use directly.

With ase (and pymatgen) the system/cell boxes always have an origin at (0,0,0), but atomman allows it to be different. This is useful for dislocations as we can put the (0,0,0) position in the center of the system with the unrelaxed dislocation core on it. The dump should be using scaled atomic positions so that the positions are automatically adjusted with the change to the box. This means that the dislocation in the ase object will not be at (0,0,0), but should now remain in the center of the system.

For the reverse, loading the ase.Atoms back into atomman shouldn't be an issue. It would just replicate the ase system's box as there is no information about where the origin should actually be.

I also updated phonopy accordingly since ase.Atoms and phonopy.Atoms are nearly identical in API. It looks like the pymatgen dump already uses the scaled positions.

Update atomman/dump/ase_Atoms/dump.py, or place in your code and call dump(system)

# coding: utf-8

# Standard Python imports
from typing import Tuple, Union

# http://www.numpy.org/
import numpy as np

# https://wiki.fysik.dtu.dk/ase/
import ase

def dump(system,
         symbols: Union[str, list, None] = None,
         return_prop: bool = False
         ) -> Union[ase.Atoms, Tuple[ase.Atoms, dict]]:
    """
    Convert an atomman.System into an ase.Atoms.

    Parameters
    ----------
    system : atomman.System
        A atomman representation of a system.
    symbols : tuple, optional
        List of the element symbols that correspond to the atom types.  If not
        given, will use system.symbols if set, otherwise no element content
        will be included.
    return_prop : bool, optional
        Indicates if the extra per-atom properties are to be returned in a
        dictionary.  Default value is False.
    Returns
    -------
    aseatoms : ase.Atoms
        An ase representation of a collection of atoms.
    prop : dict
        Dictionary containing any extra per-atom properties to include.
        Returned if return_prop is True.
    """

    # Get box/cell information
    cell = system.box.vects
    pbc = system.pbc

    # Get symbols information
    if symbols is None:
        symbols = system.symbols
    symbols = np.asarray(symbols)
    if None in symbols:
        raise ValueError('Symbols needed for all atypes')

    # Convert short symbols list to full allsymbols list
    atype = system.atoms.atype
    allsymbols = symbols[atype-1]

    # Get atomic information
    spos = system.atoms_prop('pos', scale=True)
    prop = {}
    for p in system.atoms_prop():
        if p != 'atype' and p != 'pos':
            prop[p] = system.atoms_prop(key=p)

    # Build Atoms
    aseatoms = ase.Atoms(symbols=allsymbols, scaled_positions=spos, pbc=pbc, cell=cell)

    if return_prop is True:
        return aseatoms, prop
    else:
        return aseatoms

[naveenmohandas]

Hi thank you, that works. I also had another question, would there be any issue when converting to Ase atoms object for relaxation and converting it back to Atomman system for calculating the differential displacement? Is there a standard way to keep constraints? at the moment I added the below section to the dump to add constraint in ase.Atoms.

def custom_dump(system,
         symbols: Union[str, list, None] = None,
         return_prop: bool = False
         ) -> Union[ase.Atoms, Tuple[ase.Atoms, dict]]:

    # Get box/cell information
    cell = system.box.vects
    pbc = system.pbc

    # Get symbols information
    if symbols is None:
        symbols = system.symbols
    symbols = np.asarray(symbols)
    if None in symbols:
        raise ValueError('Symbols needed for all atypes')

    # Convert short symbols list to full allsymbols list
    atype = system.atoms.atype
    allsymbols = symbols[atype-1]

    # get mask positions
    mask = get_mask(atype)

    # Get atomic information
    positions = system.atoms.pos
    spos = system.atoms_prop('pos', scale=True)
    prop = {}
    for p in system.atoms_prop():
        if p != 'atype' and p != 'pos':
            prop[p] = system.atoms_prop(key=p)

    # Build Atoms
    # using scaled position
    aseatoms = ase.Atoms(symbols=allsymbols, scaled_positions=spos, pbc=pbc, cell=cell)

    # set constraint for atoms
    from ase.constraints import FixAtoms

    constraint = FixAtoms(mask=mask)
    aseatoms.set_constraint(constraint)

    if return_prop is True:
        return aseatoms, prop
    else:
        return aseatoms

def get_mask(atype: list[int]) -> list[bool]:
    mask = []
    for i in atype:
        if i == 1:
            mask.append(False) 
        else:
            mask.append(True)
    return mask

[lmhale99]

The conversion itself should not be an issue, but when you calculate the differential displacement the relaxed system and the original base (dislocation-free) system need to have compatible atom ids.

• The Dislocation method returns both the base and unrelaxed dislocation systems.
• Converting to ase and back again will rigidly shift all atoms due to the origin issue.
• Either the base system needs to be rigidly shifted as well (convert to ase and back again), or the relaxed dislocation system needs to be rigidly shifted back before computing the differential displacement [something like system.box_set(vects=system.box.vects, origin=original_origin, scale=True) should work].
• I haven't done much ase runs, but the relaxation needs to keep the atoms in the same order and ideally not wrap atomic positions around the periodic boundary.
• Finally, when you get to the differential displacement, if the line length direction is very thin you may need to replicate along that direction using supersize to avoid atoms neighboring the same atom more than once.

With the constraints, this seems to be outside the dump behavior as it is a specific case. However, you should be able to do it fairly simply with something like

aseatoms = system.dump('ase_Atoms')

mask = (system.atoms.atype != 1).tolist()     # this works as atom properties are numpy arrays.
constraint = FixAtoms(mask=mask)
aseatoms.set_constraint(constraint)

[naveenmohandas]

Thank you! I will keep it mind when working doing the conversions.