import sisl as si
import numpy as np

Reading/writing geometries

sisl provides an interface for various electronic structure codes as well as commonly found data standards.

---

In this tutorial we will create some geometries, write them in various formats, and re-read them in.

BN = si.geom.bilayer(1.42, si.Atom("N"), si.Atom("B"))
print(BN)

Geometry{na: 4, no: 4,
 Atoms{species: 2,
  Atom{N, Z: 7, mass(au): 14.00670, maxR: -1.00000,
   Orbital{R: -1.00000, q0: 0.0}
  }: 2,
  Atom{B, Z: 5, mass(au): 10.81100, maxR: -1.00000,
   Orbital{R: -1.00000, q0: 0.0}
  }: 2,
 },
 maxR: -1.00000,
 Lattice{nsc: [3 3 1],
  origin=[0.0000, 0.0000, 0.0000],
  A=[2.1300, -1.2298, 0.0000],
  B=[2.1300, 1.2298, 0.0000],
  C=[0.0000, 0.0000, 23.3500],
  bc=[Periodic,
      Periodic,
      Unknown]
 }
}

We will now write this geometry out into a common xyz file format. This file-format is easily parseable by a large number of other codes (including being compatible with the ASE format).

BN.write("BN.xyz")

By now there should be a file called BN.xyz in the current directory.

---

Let us try and read it in again, to check it retains the same properties. There are various ways to do this: - si.io.get_sile(filename).read_geometry() - si.Geometry.read(filename) will internally do as above - si.Geometry.new(filename), this last method is the most versatile method as it can also work on Python objects

BN2 = si.Geometry.new("BN.xyz")
assert BN2 == BN

For a larger interaction with the file content, say if the file contains both geometries and real space quantities, it can be benificial to store the file handle. In sisl, files are called Sile.

xyz = si.io.get_sile("BN.xyz")

Now first read the Lattice:

print(xyz.read_lattice())

Lattice{nsc: [3 3 1],
 origin=[0.0000, 0.0000, 0.0000],
 A=[2.1300, -1.2298, 0.0000],
 B=[2.1300, 1.2298, 0.0000],
 C=[0.0000, 0.0000, 23.3500],
 bc=[Periodic,
     Periodic,
     Unknown]
}

And then read the Geometry

print(xyz.read_geometry())

Geometry{na: 4, no: 4,
 Atoms{species: 2,
  Atom{N, Z: 7, mass(au): 14.00670, maxR: -1.00000,
   Orbital{R: -1.00000, q0: 0.0}
  }: 2,
  Atom{B, Z: 5, mass(au): 10.81100, maxR: -1.00000,
   Orbital{R: -1.00000, q0: 0.0}
  }: 2,
 },
 maxR: -1.00000,
 Lattice{nsc: [3 3 1],
  origin=[0.0000, 0.0000, 0.0000],
  A=[2.1300, -1.2298, 0.0000],
  B=[2.1300, 1.2298, 0.0000],
  C=[0.0000, 0.0000, 23.3500],
  bc=[Periodic,
      Periodic,
      Unknown]
 }
}

Other file formats

There are a broad range of file formats. To automatically write out into the Siesta XV file format, simply do:

BN.write("BN.XV")

One cannot expect all file-formats to retain all information in a geometry. For instance the xyz file format does not specify how orbitals should be described. Therefore orbital information will be lost when writing to the xyz file format, see for instance here:

BN2 = si.geom.bilayer(1.42, si.Atom("N", [1, 2]), si.Atom("B", [2, 3]))
print(BN2.atoms)

Atoms{species: 2,
 Atom{N, Z: 7, mass(au): 14.00670, maxR: 2.00000,
  Orbital{R: 1.00000, q0: 0.0},
  Orbital{R: 2.00000, q0: 0.0}
 }: 2,
 Atom{B, Z: 5, mass(au): 10.81100, maxR: 3.00000,
  Orbital{R: 2.00000, q0: 0.0},
  Orbital{R: 3.00000, q0: 0.0}
 }: 2,
}

Note how the two atoms has multiple orbitals, with different orbital ranges.

BN2.write("BN2.xyz")
print(si.Geometry.new("BN2.xyz").atoms)

Atoms{species: 2,
 Atom{N, Z: 7, mass(au): 14.00670, maxR: -1.00000,
  Orbital{R: -1.00000, q0: 0.0}
 }: 2,
 Atom{B, Z: 5, mass(au): 10.81100, maxR: -1.00000,
  Orbital{R: -1.00000, q0: 0.0}
 }: 2,
}

Here we extracted only the atoms object to show the difference there.

Selecting origin of the output

sisl implements a variety of objects that interacts with the stdout of codes. For instance:

# Regular siesta out
siesta RUN.fdf > RUN.out
# VASP out
vasp > RUN.out
...

In general the extensions are not well-defined and there is a high probability of overlapping extensions with different codes. To make it simpler for the user to use the correct object for a file, one can specify a name of the code origin:

siesta_out = si.get_sile("RUN.out{siesta}")
vasp_out = si.get_sile("RUN.out{vasp}")