Skip to content

atomlib.defect

A collection of functions for inserting dislocations into structures.

CutType module-attribute

CutType: TypeAlias = Literal['shift', 'add', 'rm']

Cut plane to use when creating a (non-screw) dislocation.

ellip_pi

ellip_pi(
    n: NDArray[float64], m: NDArray[float64]
) -> NDArray[float64]

Complete elliptic integral of the third kind, \(\Pi(n | m)\).

Follows the definition of Wolfram Mathworld.

Source code in atomlib/defect.py
def ellip_pi(n: NDArray[numpy.float64], m: NDArray[numpy.float64]) -> NDArray[numpy.float64]:
    """
    Complete elliptic integral of the third kind, $\\Pi(n | m)$.

    Follows the definition of [Wolfram Mathworld][wolfram_ellip_pi].

    [wolfram_ellip_pi]: https://mathworld.wolfram.com/EllipticIntegraloftheThirdKind.html
    """
    from scipy.special import elliprf, elliprj  # type: ignore

    y = 1 - m
    assert numpy.all(y > 0)

    rf = elliprf(0, y, 1)
    rj = elliprj(0, y, 1, 1 - n)
    return rf + rj * n / 3

stacking_fault

stacking_fault(
    atoms: HasAtomsT,
    pos: VecLike,
    shift: VecLike,
    plane: VecLike,
) -> HasAtomsT

Add a stacking fault to the structure.

The fault plane will pass through the position pos, with normal plane. Atoms above plane will be shifted by the vector shift.

If there is a component of shift parallel to plane, applying the shift will create (or remove) volume. In this case, atoms are added (or removed) accordingly. If shift \(\cdot\) plane is positive, atoms are added.

In general, adding a stacking fault will not preserve a cell's periodicity. Please verify this for your usecase.

PARAMETER DESCRIPTION
atoms

Structure to add fault to

TYPE: HasAtomsT

pos

Position on fault plane

TYPE: VecLike

shift

Vector to shift by

TYPE: VecLike

plane

Normal to fault plane

TYPE: VecLike

RETURNS DESCRIPTION
HasAtomsT

Structure with a stacking fault added

Source code in atomlib/defect.py
def stacking_fault(atoms: HasAtomsT, pos: VecLike, shift: VecLike, plane: VecLike) -> HasAtomsT:
    """
    Add a stacking fault to the structure.

    The fault plane will pass through the position `pos`, with normal `plane`.
    Atoms above `plane` will be shifted by the vector `shift`.

    If there is a component of `shift` parallel to `plane`, applying the shift
    will create (or remove) volume. In this case, atoms are added (or removed)
    accordingly. If `shift` $\\cdot$ `plane` is positive, atoms are added.

    In general, adding a stacking fault will not preserve a cell's periodicity.
    Please verify this for your usecase.

    Args:
      atoms: Structure to add fault to
      pos: Position on fault plane
      shift: Vector to shift by
      plane: Normal to fault plane

    Returns:
      Structure with a stacking fault added
    """
    pos = to_vec3(pos)
    shift = to_vec3(shift)
    plane = to_vec3(plane)
    plane /= numpy.linalg.norm(plane)

    coords = atoms.coords(frame='local')  # type: ignore

    perp_dist = numpy.dot(coords - pos, plane)
    atoms = atoms.with_columns(perp_dist=polars.Series(perp_dist)) \
        .with_columns(branch=polars.col('perp_dist') > 0.)

    d = numpy.dot(shift, plane)
    if -d > 1e-8:
        logging.info("Removing atoms.")
        old_len = len(atoms)
        atoms = atoms.filter(
            (polars.col('perp_dist') < 0) | (polars.col('perp_dist') >= -d)
        )
        logging.info(f"Removed {old_len - len(atoms)} atoms")
    if d > 1e-8:
        logging.info("Duplicating atoms.")
        duplicate = atoms.filter(
            (polars.col('perp_dist') >= -d) & (polars.col('perp_dist') < 0)
        ).with_columns(~polars.col('branch'))
        logging.info(f"Duplicated {len(duplicate)} atoms")

        atoms = atoms.with_atoms(Atoms.concat((atoms, duplicate)))

    coords = atoms.coords(frame='local')
    branch = atoms['branch'].to_numpy()
    atoms = atoms.with_coords(coords + shift * branch[:, None], frame='local')
    return atoms.with_atoms(Atoms(atoms.drop('perp_dist', 'branch')))

disloc_edge

disloc_edge(
    atoms: HasAtomsT,
    center: VecLike,
    b: VecLike,
    t: VecLike,
    cut: Union[CutType, VecLike] = "shift",
    *,
    poisson: float = 0.25
) -> HasAtomsT

Add a Volterra edge dislocation to the structure.

The dislocation will pass through center, with line vector t and Burgers vector b. t will be modified such that it is perpendicular to b.

The cut parameter defines the cut (discontinuity) plane used to displace the atoms. By default, cut is 'shift', which defines the cut plane to contain b, ensuring no atoms need to be added or removed. Instead, 'add' or 'rm' may be specified, which defines the cut plane as containing \(\mathbf{b} \times \mathbf{t}\). In this mode, atoms will be added or removed in the plane of b to create the dislocation. Alternatively, a vector cut may be supplied. The cut plane will contain this vector.

In the coordinate system where t is along \(z\), and b is along \(+x\), the displacement due to the edge dislocation can be calculated as follows:

\[\begin{align} u_x &= \frac{b}{2\pi} \left( \arctan(x, y) + \frac{x y}{2(x^2 + y^2)(1-\nu)} \right) \\ u_y &= -\frac{b}{4(1-\nu)} \left( (1-2\nu) \ln(x^2 + y^2) + \frac{x^2 - y^2}{x^2 + y^2} \right) \end{align}\]

Where \(x\) and \(y\) are distances from the dislocation center. This creates a discontinuity along the \(-x\) axis. This coordinate system is rotated to support branches along an arbitrary axis.

The dislocation is defined by the FS/RH convention: Performing one loop of positive sense relative to the tangent vector displaces the real crystal one b relative to a reference crystal.

PARAMETER DESCRIPTION
atoms

Structure to add edge dislocation to

TYPE: HasAtomsT

center

Position on dislocation line

TYPE: VecLike

b

Burgers vector of dislocation (FS/RH convention)

TYPE: VecLike

t

Tangent vector of dislocation

TYPE: VecLike

cut

Cut plane to create dislocation on

TYPE: Union[CutType, VecLike] DEFAULT: 'shift'

poisson

Poisson ratio of material

TYPE: float DEFAULT: 0.25

RETURNS DESCRIPTION
HasAtomsT

Structure with an edge dislocation added

Source code in atomlib/defect.py
def disloc_edge(atoms: HasAtomsT, center: VecLike, b: VecLike, t: VecLike, cut: t.Union[CutType, VecLike] = 'shift',
                *, poisson: float = 0.25) -> HasAtomsT:
    r"""
    Add a Volterra edge dislocation to the structure.

    The dislocation will pass through `center`, with line vector `t` and Burgers vector `b`.
    `t` will be modified such that it is perpendicular to `b`.

    The `cut` parameter defines the cut (discontinuity) plane used to displace the atoms.
    By default, `cut` is `'shift'`, which defines the cut plane to contain `b`, ensuring no atoms
    need to be added or removed.
    Instead, `'add'` or `'rm'` may be specified, which defines the cut plane as containing
    $\mathbf{b} \times \mathbf{t}$. In this mode, atoms will be added or removed
    in the plane of `b` to create the dislocation. Alternatively, a vector `cut`
    may be supplied. The cut plane will contain this vector.

    In the coordinate system where `t` is along $z$, and `b` is along $+x$, the displacement
    due to the edge dislocation can be calculated as follows:

    $$\begin{align}
       u_x &= \frac{b}{2\pi} \left( \arctan(x, y) + \frac{x y}{2(x^2 + y^2)(1-\nu)} \right) \\
       u_y &= -\frac{b}{4(1-\nu)} \left( (1-2\nu) \ln(x^2 + y^2) + \frac{x^2 - y^2}{x^2 + y^2} \right)
    \end{align}$$

    Where $x$ and $y$ are distances from the dislocation center. This creates a discontinuity along
    the $-x$ axis. This coordinate system is rotated to support branches along an arbitrary axis.

    The dislocation is defined by the FS/RH convention: Performing one loop of positive sense
    relative to the tangent vector displaces the real crystal one `b` relative to a reference
    crystal.

    Args:
      atoms: Structure to add edge dislocation to
      center: Position on dislocation line
      b: Burgers vector of dislocation (FS/RH convention)
      t: Tangent vector of dislocation
      cut: Cut plane to create dislocation on
      poisson: Poisson ratio of material

    Returns:
      Structure with an edge dislocation added
    """

    center = to_vec3(center)
    b_vec = to_vec3(b)
    #b_mag = norm(b_vec)

    # get component of t perpendicular to b, normalize
    t = to_vec3(t)
    t = perp(t, b)
    if norm(t) < 1e-10:
        raise ValueError("`b` and `t` must be different.")
    t /= norm(t)

    if isinstance(cut, str):
        cut = cast(CutType, cut.lower())
        if cut == 'shift':
            plane_v = b_vec.copy()
        elif cut == 'add':
            # FS/RH convention: t x b points to extra half plane
            plane_v = numpy.cross(t, b_vec)
        elif cut == 'rm':
            plane_v = -numpy.cross(t, b_vec)
        else:
            raise ValueError(f"Unknown cut plane type `{cut}`. Expected 'shift', 'add', 'rm', or a vector.")
        plane_v /= norm(plane_v)
    else:
        plane_v = to_vec3(cut)
        plane_v = plane_v / norm(plane_v)
        if numpy.linalg.norm(numpy.cross(plane_v, t)) < 1e-10:
            raise ValueError('`cut` and `t` must be different.')

    # translate center to 0., and align t to [0, 0, 1], plane to +y
    transform = AffineTransform3D.translate(center).inverse().compose(
        LinearTransform3D.align_to(t, [0., 0., 1.], plane_v, [-1., 0., 0.])
    )
    frame = atoms.get_atoms('local').transform(transform)
    b_vec = transform.transform_vec(b_vec)

    d = numpy.dot(b_vec, [0., 1., 0.])
    if -d > 1e-8:
        logging.info("Removing atoms.")
        old_len = len(frame)
        frame = frame.filter(~(
            (polars.col('coords').arr.get(0) < 0)
            & (polars.col('coords').arr.get(1) >= d/2.)
            & (polars.col('coords').arr.get(1) <= -d/2.)
        ))
        logging.info(f"Removed {old_len - len(frame)} atoms")
    if d > 1e-8:
        logging.info("Duplicating atoms.")
        duplicate = frame.filter(
            (polars.col('coords').arr.get(0) < 0)
            & (polars.col('coords').arr.get(1) >= -d/2.)
            & (polars.col('coords').arr.get(1) <= d/2.)
        )
        logging.info(f"Duplicated {len(duplicate)} atoms")

        frame = Atoms.concat((frame, duplicate))
        #atoms = atoms._replace_atoms(frame)
        branch = numpy.ones(len(frame), dtype=float)
        if len(duplicate) > 0:
            branch[-len(duplicate):] *= -1  # flip branch of duplicated atoms
    else:
        branch = numpy.ones(len(frame), dtype=float)

    pts = frame.coords()

    x, y, _ = split_arr(pts, axis=-1)
    r2 = x**2 + y**2

    # displacement parallel and perpendicular to b
    d_para = branch * numpy.arctan2(y, x) + x*y/(2*(1-poisson)*r2)
    d_perp = -(1-2*poisson)/(4*(1-poisson)) * numpy.log(r2) + (y**2 - x**2)/(4*(1-poisson)*r2)

    disps = numpy.stack([
        d_para * b_vec[0] + d_perp * b_vec[1],
        d_perp * b_vec[0] + d_para * b_vec[1],
        numpy.zeros_like(x)
    ], axis=-1) / (2*numpy.pi)

    return atoms.with_atoms(frame.with_coords(pts + disps).transform(transform.inverse()), 'local')

disloc_screw

disloc_screw(
    atoms: HasAtomsT,
    center: VecLike,
    b: VecLike,
    cut: Optional[VecLike] = None,
    sign: bool = True,
) -> HasAtomsT

Add a Volterra screw dislocation to the structure.

The dislocation will pass through center, with Burgers vector b.

The cut parameter defines the cut (discontinuity) plane used to displace the atoms. By default, cut is chosen automtically, but it may also be specified as a vector which points from the dislocation core towards the cut plane (not normal to the cut plane!)

The screw dislocation in an isotropic medium has a particularily simple form, given by:

\[ \mathbf{u} = \frac{\mathbf{b}}{2\pi} \arctan(x, y) \]

for a dislocation along \(+z\) with cut plane along \(-x\). To support arbitrary cut planes, \(x\) and \(y\) are replaced by the components of \(r\) parallel and perpendicular to the cut plane, evaluated in the plane of b.

The dislocation is defined by the FS/RH convention: Performing one loop of positive sense relative to the tangent vector displaces the real crystal one b relative to a reference crystal.

PARAMETER DESCRIPTION
atoms

Structure to add screw dislocation to

TYPE: HasAtomsT

center

Position on dislocation line

TYPE: VecLike

b

Burgers vector of dislocation (FS/RH convention)

TYPE: VecLike

cut

Cut plane to create dislocation on

TYPE: Optional[VecLike] DEFAULT: None

sign

Sign of screw dislocation (True means positive)

TYPE: bool DEFAULT: True

RETURNS DESCRIPTION
HasAtomsT

Structure with a screw dislocation added

Source code in atomlib/defect.py
def disloc_screw(atoms: HasAtomsT, center: VecLike, b: VecLike, cut: t.Optional[VecLike] = None,
                 sign: bool = True) -> HasAtomsT:
    r"""
    Add a Volterra screw dislocation to the structure.

    The dislocation will pass through `center`, with Burgers vector `b`.

    The `cut` parameter defines the cut (discontinuity) plane used to displace the atoms.
    By default, `cut` is chosen automtically, but it may also be specified as a vector
    which points from the dislocation core towards the cut plane (not normal to the cut plane!)

    The screw dislocation in an isotropic medium has a particularily simple form, given by:

    $$
       \mathbf{u} = \frac{\mathbf{b}}{2\pi} \arctan(x, y)
    $$

    for a dislocation along $+z$ with cut plane along $-x$. To support arbitrary cut planes,
    $x$ and $y$ are replaced by the components of $r$ parallel and perpendicular to the cut plane,
    evaluated in the plane of `b`.

    The dislocation is defined by the FS/RH convention: Performing one loop of positive sense
    relative to the tangent vector displaces the real crystal one `b` relative to a reference
    crystal.

    Args:
      atoms: Structure to add screw dislocation to
      center: Position on dislocation line
      b: Burgers vector of dislocation (FS/RH convention)
      cut: Cut plane to create dislocation on
      sign: Sign of screw dislocation (True means positive)

    Returns:
      Structure with a screw dislocation added
    """

    center = to_vec3(center)
    b_vec = to_vec3(b)
    t = b_vec / float(numpy.linalg.norm(b_vec))
    t = -t if not sign else t
    if cut is None:
        if numpy.linalg.norm(numpy.cross(t, [1., 1., 1.])) < numpy.pi/4:
            # near 111, choose x as cut plane direction
            cut = to_vec3([1., 0., 0.])
        else:
            # otherwise find plane by rotating around 111
            cut = cast(NDArray[numpy.float64], LinearTransform3D.rotate([1., 1., 1.], 2*numpy.pi/3).transform(t))
    else:
        cut = to_vec3(cut)
        cut /= norm(cut)
        if numpy.allclose(cut, t, atol=1e-2):
            raise ValueError("`t` and `cut` must be different.")

    print(f"Cut plane direction: {cut}")

    frame = atoms.get_atoms('local')
    pts = frame.coords() - center

    # components perpendicular to t
    cut_perp = -perp(cut, t)
    pts_perp = perp(pts, t)

    # signed angle around dislocation
    theta = numpy.arctan2(dot(t, numpy.cross(cut_perp, pts_perp)), dot(cut_perp, pts_perp))
    # FS/RH convention
    disp = b_vec * (theta / (2*numpy.pi))

    return atoms.with_atoms(frame.with_coords(pts + center + disp), 'local')

disloc_loop_z

disloc_loop_z(
    atoms: HasAtomsT,
    center: VecLike,
    b: VecLike,
    loop_r: float,
    *,
    poisson: float = 0.25
) -> HasAtomsT

Add a circular dislocation loop to the structure, assuming an elastically isotropic material.

The loop will have radius loop_r, be centered at center, and oriented along the z-axis.

The dislocation's sign is defined such that travelling upwards through the loop results in a displacement of b. poisson is the material's poisson ratio, which affects the shape of dislocations with an edge component.

Adding the loop creates (or removes) a volume of \(\mathbf{b} \cdot \mathbf{n}A\), where \(\mathbf{n}A\) is the loop's normal times its area. Consequently, this function adds or remove atoms to effect this volume change.

PARAMETER DESCRIPTION
atoms

Structure to add dislocation loop to

TYPE: HasAtomsT

center

Center of dislocation loop

TYPE: VecLike

b

Burgers vector of dislocation (defined so travelling upwards leads to a plastic displacment of b).

TYPE: VecLike

loop_r

Radius of dislocation loop

TYPE: float

poisson

Poisson ratio of material

TYPE: float DEFAULT: 0.25

RETURNS DESCRIPTION
HasAtomsT

Structure with a dislocation loop added

Source code in atomlib/defect.py
def disloc_loop_z(atoms: HasAtomsT, center: VecLike, b: VecLike,
                  loop_r: float, *, poisson: float = 0.25) -> HasAtomsT:
    r"""
    Add a circular dislocation loop to the structure, assuming an elastically isotropic material.

    The loop will have radius `loop_r`, be centered at `center`, and oriented along the z-axis.

    The dislocation's sign is defined such that travelling upwards through the loop results in a displacement of `b`.
    `poisson` is the material's poisson ratio, which affects the shape of dislocations with an edge component.

    Adding the loop creates (or removes) a volume of $\mathbf{b} \cdot \mathbf{n}A$, where $\mathbf{n}A$ is the loop's
    normal times its area. Consequently, this function adds or remove atoms to effect this volume change.

    Args:
      atoms: Structure to add dislocation loop to
      center: Center of dislocation loop
      b: Burgers vector of dislocation (defined so travelling upwards leads to a plastic displacment of `b`).
      loop_r: Radius of dislocation loop
      poisson: Poisson ratio of material

    Returns:
      Structure with a dislocation loop added
    """

    center = to_vec3(center)
    b_vec = to_vec3(b)

    atoms = atoms.transform_atoms(AffineTransform3D.translate(center).inverse())
    frame = atoms.get_atoms('local')
    branch = None

    d = numpy.dot(b_vec, [0, 0, 1])
    if -d > 1e-8:
        logging.info("Non-conservative dislocation. Removing atoms.")
        frame = frame.filter(~(
            (frame.x()**2 + frame.y()**2 < loop_r**2)
            & frame.z().is_between(d/2., -d/2., closed='both')
        ))

    if d > 1e-8:
        logging.info("Non-conservative dislocation. Duplicating atoms.")
        duplicate = frame.filter(
            (frame.x()**2 + frame.y()**2 < loop_r**2)
            & frame.z().is_between(-d/2., d/2., closed='both')
        )
        logging.info(f"Adding {len(duplicate)} atoms")

        frame = Atoms.concat((frame, duplicate))
        #atoms = atoms._replace_atoms(frame)
        branch = numpy.sign(frame['z'].to_numpy())
        if len(duplicate) > 0:
            branch[-len(duplicate):] *= -1  # flip branch of duplicated atoms

    pts = frame.coords()
    disps = _loop_disp_z(pts, b_vec, loop_r, poisson=poisson, branch=branch)

    return atoms.with_atoms(frame.with_coords(pts + disps + center), 'local')

disloc_square_z

disloc_square_z(
    atoms: HasAtomsT,
    center: VecLike,
    b: VecLike,
    loop_r: float,
    *,
    poisson: float = 0.25
) -> HasAtomsT

Add a square dislocation loop to the structure, assuming an elastically isotropic material.

The dislocation's sign is defined such that traveling upwards through the loop results in a displacement of b. poisson is the material's poisson ratio, which affects the shape of dislocations with an edge component.

Adding the loop creates (or removes) a volume of \(\mathbf{b} \cdot \mathbf{n}A\), where \(\mathbf{n}A\) is the loop's normal times its area. Consequently, this function adds or remove atoms to effect this volume change.

PARAMETER DESCRIPTION
atoms

Structure to add dislocation loop to

TYPE: HasAtomsT

center

Center of dislocation loop

TYPE: VecLike

b

Burgers vector of dislocation (defined so travelling upwards leads to a plastic displacment of b).

TYPE: VecLike

loop_r

Radius (side length/2) of dislocation loop

TYPE: float

poisson

Poisson ratio of material

TYPE: float DEFAULT: 0.25

RETURNS DESCRIPTION
HasAtomsT

Structure with a dislocation loop added

Source code in atomlib/defect.py
def disloc_square_z(atoms: HasAtomsT, center: VecLike, b: VecLike,
                    loop_r: float, *, poisson: float = 0.25) -> HasAtomsT:
    r"""
    Add a square dislocation loop to the structure, assuming an elastically isotropic material.

    The dislocation's sign is defined such that traveling upwards through the loop results in a displacement of `b`.
    `poisson` is the material's poisson ratio, which affects the shape of dislocations with an edge component.

    Adding the loop creates (or removes) a volume of $\mathbf{b} \cdot \mathbf{n}A$, where $\mathbf{n}A$ is the loop's
    normal times its area. Consequently, this function adds or remove atoms to effect this volume change.

    Args:
      atoms: Structure to add dislocation loop to
      center: Center of dislocation loop
      b: Burgers vector of dislocation (defined so travelling upwards leads to a plastic displacment of `b`).
      loop_r: Radius (side length/2) of dislocation loop
      poisson: Poisson ratio of material

    Returns:
      Structure with a dislocation loop added
    """
    poly = loop_r * numpy.array([(1, 1), (-1, 1), (-1, -1), (1, -1)])
    return disloc_poly_z(atoms, b, poly, center, poisson=poisson)

disloc_poly_z

disloc_poly_z(
    atoms: HasAtomsT,
    b: VecLike,
    poly: ArrayLike,
    center: Optional[VecLike] = None,
    *,
    poisson: float = 0.25
) -> HasAtomsT

Add a dislocation loop defined by the polygon poly, assuming an elastically isotropic material.

poly should be a 2d polygon (shape (N, 2)). It will be placed at center, in the plane z=center[2]. For CCW winding order, traveling upwards through the plane of the loop results in a displacement of b. poisson is the material's poisson ratio, which affects the shape of dislocations with an edge component.

Adding the loop creates (or removes) a volume of \(\mathbf{b} \cdot \mathbf{n}A\), where \(\mathbf{n}A\) is the loop's normal times its area. Consequently, this function adds or remove atoms to effect this volume change.

Follows the solution in Hirth, J. P. & Lothe, J. (1982). Theory of Dislocations. ISBN 0-89464-617-6

PARAMETER DESCRIPTION
atoms

Structure to add dislocation loop to

TYPE: HasAtomsT

b

Burgers vector of dislocation (defined so travelling upwards leads to a plastic displacment of b).

TYPE: VecLike

poly

2D polygon defining dislocation line. It is placed in the plane z=center[2].

TYPE: ArrayLike

center

Center of dislocation loop (defaults to zero, applying no displacement to the gven polygon).

TYPE: Optional[VecLike] DEFAULT: None

poisson

Poisson ratio of material

TYPE: float DEFAULT: 0.25

RETURNS DESCRIPTION
HasAtomsT

Structure with a dislocation loop added

Source code in atomlib/defect.py
def disloc_poly_z(atoms: HasAtomsT, b: VecLike, poly: ArrayLike, center: t.Optional[VecLike] = None,
                  *, poisson: float = 0.25) -> HasAtomsT:
    r"""
    Add a dislocation loop defined by the polygon `poly`, assuming an elastically isotropic material.

    `poly` should be a 2d polygon (shape `(N, 2)`). It will be placed at `center`, in the plane `z=center[2]`.
    For CCW winding order, traveling upwards through the plane of the loop results in a displacement of `b`.
    `poisson` is the material's poisson ratio, which affects the shape of dislocations with an edge component.

    Adding the loop creates (or removes) a volume of $\mathbf{b} \cdot \mathbf{n}A$, where $\mathbf{n}A$ is the loop's
    normal times its area. Consequently, this function adds or remove atoms to effect this volume change.

    Follows the solution in [Hirth, J. P. & Lothe, J. (1982). Theory of Dislocations. ISBN 0-89464-617-6
    ](https://www.google.com/books/edition/Theory_of_Dislocations/TAjwAAAAMAAJ)

    Args:
      atoms: Structure to add dislocation loop to
      b: Burgers vector of dislocation (defined so travelling upwards leads to a plastic displacment of `b`).
      poly: 2D polygon defining dislocation line. It is placed in the plane `z=center[2]`.
      center: Center of dislocation loop (defaults to zero, applying no displacement to the gven polygon).
      poisson: Poisson ratio of material

    Returns:
      Structure with a dislocation loop added
    """
    center = to_vec3(center if center is not None else [0., 0., 0.])
    b_vec = to_vec3(b)

    poly = numpy.atleast_2d(poly)
    if poly.ndim != 2 or poly.shape[-1] != 2:
        raise ValueError(f"Expected a 2d polygon. Instead got shape {poly.shape}")

    frame = atoms.get_atoms('local')
    coords: NDArray[numpy.float64] = frame.coords() - center

    branch = None
    d = numpy.dot(b_vec, [0, 0, 1])
    if abs(d) > 1e-8:
        logging.info("Non-conservative dislocation.")
        windings = polygon_winding(poly, coords[..., :2])

        z = coords[..., 2]
        remove = (z >= windings * d/2.) & (z <= -windings * d/2.)
        duplicate = (z >= -windings * d/2.) & (z <= windings * d/2.)

        n_remove = numpy.sum(remove, dtype=int)
        if n_remove:
            logging.info(f"Removing {n_remove} atoms")
            frame = frame.filter(polars.Series(~remove))
            duplicate = duplicate[~remove]

        n_duplicate = numpy.sum(duplicate, dtype=int)
        if n_duplicate:
            logging.info(f"Duplicating {n_duplicate} atoms")
            frame = Atoms.concat((frame, frame.filter(polars.Series(duplicate))))

            branch = numpy.ones(len(frame))
            branch[-n_duplicate:] = -1  # flip branch of duplicated atoms

        coords = frame.coords() - center

    disp = _poly_disp_z(coords, b_vec, poly, poisson=poisson, branch=branch)

    return atoms.with_atoms(frame.with_coords(coords + disp + center), 'local')