# Creating structures in pyiron¶

This section gives a brief introduction about some of the tools available in pyiron to construct atomic structures.

For the sake of compatibility, our structure class is written to be compatible with the popular Atomistic Simulation Environment package (ASE). This makes it possible to use routines from ASE to help set-up structures.

Furthermore, pyiron uses the NGLview package to visualize the structures and trajectories interactively in 3D using NGLview-widgets.

As preparation for the following discussion we import a few python libraries

In [1]:
import numpy as np
%matplotlib inline
import matplotlib.pylab as plt


and create a pyiron project named 'structures':

In [2]:
from pyiron.project import Project
pr = Project(path='structures')


## Bulk crystals¶

In this section we discuss various possibilities to create bulk crystal structures.

### Using create_structure()¶

The simplest way to generate simple crystal structures is using the inbuilt create_structure() function specifying the element symbol, Bravais basis and the lattice constant(s)

Note: The output gives a cubic cell rather than the smallest non-orthogonal unit cell.

In [3]:
structure = pr.create_structure('Al',
bravais_basis='fcc',
lattice_constant=4.05)


To plot the structure interactively in 3D simply use:

In [4]:
structure.plot3d()


### Using create_ase_bulk()¶

Another convenient way to set up structures is using the create_ase_bulk() function which is built on top of the ASE build package for bulk crystals. This function returns an object which is of the pyiron structure object type.

Example: fcc bulk aluminum in a cubic cell

In [5]:
structure = pr.create_ase_bulk('Al', cubic=True)
structure.plot3d()


Example: wurtzite GaN in a 3x3x3 repeated orthorhombic cell.

Note:

• In contrast to new_structure = structure.repeat() which creates a new object, set_repeat() modifies the existing structure object.
• Setting spacefill=False in the plot3d() method changes the atomic structure style to "ball and stick".
In [6]:
structure = pr.create_ase_bulk('AlN',
crystalstructure='wurtzite',
a=3.5, orthorhombic=True)
structure.set_repeat([3,3,3])
structure.plot3d(spacefill=False)


## Creating surfaces (using ASE)¶

Surfaces can be created using the create_surface() function which is also built on top of the ASE build package for surfaces

Example: Creating a 3x4 fcc Al(111) surface with 4 layers and a vacuum of 10 Ångström

In [7]:
Al_111 = pr.create_surface("Al", surface_type="fcc111",
size=(3, 4, 4), vacuum=10, orthogonal=True)
Al_111.plot3d()


## Creating structures without importing the project class¶

In all the examples shown above, the structures are create from the pyiron Project object. It is also possible to do this without importing/initializing this object. For this the appropriate imports must be made.

In [8]:
from pyiron import create_ase_bulk, create_surface

In [9]:
structure = create_ase_bulk('AlN',
crystalstructure='wurtzite',
a=3.5, orthorhombic=True)
structure.set_repeat([3,3,3])
structure.plot3d(spacefill=False)

In [10]:
Al_111 = create_surface("Al", surface_type="fcc111",
size=(3, 4, 4), vacuum=10, orthogonal=True)
Al_111.plot3d()


### Using the ASE spacegroup class¶

In [11]:
from ase.spacegroup import crystal
from pyiron import ase_to_pyiron

a = 9.04
skutterudite = crystal(('Co', 'Sb'),
basis=[(0.25, 0.25, 0.25), (0.0, 0.335, 0.158)],
spacegroup=204,
cellpar=[a, a, a, 90, 90, 90])
skutterudite = ase_to_pyiron(skutterudite)

In [12]:
skutterudite.plot3d()


## Accessing the properties of the structure object¶

Using the bulk aluminum fcc example from before the structure object can be created by

In [13]:
structure = pr.create_ase_bulk('Al', cubic=True)


A summary of the information about the structure is given by using

In [14]:
print(structure)

Al: [0. 0. 0.]
Al: [0.    2.025 2.025]
Al: [2.025 0.    2.025]
Al: [2.025 2.025 0.   ]
pbc: [ True  True  True]
cell:
[[4.05 0.   0.  ]
[0.   4.05 0.  ]
[0.   0.   4.05]]



The cell vectors of the structure object can be accessed and edited through

In [15]:
structure.cell

Out[15]:
array([[4.05, 0.  , 0.  ],
[0.  , 4.05, 0.  ],
[0.  , 0.  , 4.05]])

The positions of the atoms in the structure object can be accessed and edited through

In [16]:
structure.positions

Out[16]:
array([[0.   , 0.   , 0.   ],
[0.   , 2.025, 2.025],
[2.025, 0.   , 2.025],
[2.025, 2.025, 0.   ]])

## Point defects¶

### Creating a single vacancy¶

We start by setting up a 4x4x4 supercell

In [17]:
structure = pr.create_ase_bulk('Al', cubic=True)
structure.set_repeat([4,4,4])


To create the vacancy at position index "0" simply use:

In [18]:
del structure[0]


To plot the structure that now contains a vacancy run:

In [19]:
structure.plot3d()


### Creating multiple vacancies¶

In [20]:
# First create a 4x4x4 supercell
structure = pr.create_ase_bulk('Al', cubic=True)
structure.set_repeat([4,4,4])
print('Number of atoms in the repeat unit: ',structure.get_number_of_atoms())

Number of atoms in the repeat unit:  256


The del command works for passing a list of indices to the structure object. For example, a random set of n$_{\text{vac}}$ vacancies can be created by using

In [21]:
# Generate a list of indices for the vacancies
n_vac = 24
vac_ind_lst = np.random.permutation(len(structure))[:n_vac]

# Remove atoms according to the "vac_ind_lst"
del structure[vac_ind_lst]

In [22]:
# Visualize the structure
print('Number of atoms in the repeat unit: ',structure.get_number_of_atoms())
structure.plot3d()

Number of atoms in the repeat unit:  232


### Random substitutial alloys¶

In [23]:
# Create a 4x4x4 supercell
structure = pr.create_ase_bulk('Al', cubic=True)
structure.set_repeat([4,4,4])


Substitutional atoms can be defined by changing the atomic species accessed through its position index.

Here, we set $n_{\text{sub}}$ magnesium substitutional atoms at random positions

In [24]:
n_sub = 24
structure[np.random.permutation(len(structure))[:n_sub]] = 'Mg'

In [25]:
# Visualize the structure and print some additional information about the structure
print('Number of atoms in the repeat unit: ',structure.get_number_of_atoms())
print('Chemical formula: ',structure.get_chemical_formula())
structure.plot3d()

Number of atoms in the repeat unit:  256
Chemical formula:  Al232Mg24


## Explicit definition of the structure¶

You can also set-up structures through the explicit input of the cell parameters and positions

In [26]:
cell = 10.0 * np.eye(3) # Specifying the cell dimensions
positions = [[0.25, 0.25, 0.25], [0.75, 0.75, 0.75]]
elements = ['O', 'O']

# Now use the Atoms class to create the instance.
O_dimer = pr.create_atoms(elements=elements, scaled_positions=positions, cell=cell)

O_dimer.plot3d()


## Importing from cif/other file formats¶

Parsers from ASE can be used to import structures from other formats. In this example, we will download and import a Nepheline structure from the Crystallography Open Database (COD)

In [27]:
# The COD structures can be accessed through their unique COD identifier
cod = 1008753
filename = '{}.cif'.format(cod)
url = 'http://www.crystallography.net/cod/{}'.format(filename)

In [28]:
# Download and save the structure file locally
import urllib
urllib.request.urlretrieve(url=url, filename='strucs.'+filename);

In [29]:
# Using ase parsers to read the structure and then convert to a pyiron instance
import ase
from pyiron import ase_to_pyiron

format='cif'))
structure.info["cod"] = cod

/home/surendralal/miniconda3/envs/pyiron_workshop/lib/python3.7/site-packages/ase/io/cif.py:375: UserWarning: crystal system 'hexagonal' is not interpreted for space group Spacegroup(173, setting=1). This may result in wrong setting!
setting_name, spacegroup))

In [30]:
structure.plot3d()

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
~/science/phd/dev/pyiron/pyiron/atomistics/structure/atoms.py in plot3d(self, show_cell, show_axes, camera, spacefill, particle_size, select_atoms, background, color_scheme, colors, scalar_field, scalar_start, scalar_end, scalar_cmap, vector_field, vector_color, magnetic_moments, custom_array, custom_3darray, view_plane, distance_from_camera)
1326         try:  # If the graphical packages are not available, the GUI will not work.
-> 1327             import nglview
1328         except ImportError:

ModuleNotFoundError: No module named 'nglview'

During handling of the above exception, another exception occurred:

ImportError                               Traceback (most recent call last)
<ipython-input-30-184096eccb97> in <module>
----> 1 structure.plot3d()

~/science/phd/dev/pyiron/pyiron/atomistics/structure/atoms.py in plot3d(self, show_cell, show_axes, camera, spacefill, particle_size, select_atoms, background, color_scheme, colors, scalar_field, scalar_start, scalar_end, scalar_cmap, vector_field, vector_color, magnetic_moments, custom_array, custom_3darray, view_plane, distance_from_camera)
1327             import nglview
1328         except ImportError:
-> 1329             raise ImportError(
1330                 "The package nglview needs to be installed for the plot3d() function!"
1331             )

ImportError: The package nglview needs to be installed for the plot3d() function!

Structures can be stored indepently from jobs in HDF5 by using the special StructureContainer job. To save to disk, call run().

In [32]:
container = pr.create_job(pr.job_type.StructureContainer, "nepheline")
container.structure = structure
container.run()

The job nepheline was saved and received the ID: 47


It's also possible to store multiple structures in one container and to store directly from a job. Let's use this here to store the equilibrated structures at finite temperatures.

In [46]:
al_container = pr.create_job(pr.job_type.StructureContainer, "al_temp", delete_existing_job=True)
for T in (400, 600, 800):
j = pr.create_job(pr.job_type.Lammps, "T_{}".format(T))
j.structure = pr.create_ase_bulk("Al", cubic = True)
j.potential = j.list_potentials()[0]
j.calc_md(temperature=T, n_ionic_steps=1000, pressure=0)
j.run()
structure = j.get_structure(-1)
structure.info["T"] = T
structure.info["P"] = 0
al_container.append(structure)

al_container.run()

This group does not exist in the HDF5 file al_temp

2020-08-22 14:05:03,496 - pyiron_log - WARNING - The job T_400 is being loaded instead of running. To re-run use the argument 'delete_existing_job=True in create_job'
2020-08-22 14:05:04,282 - pyiron_log - WARNING - The job T_600 is being loaded instead of running. To re-run use the argument 'delete_existing_job=True in create_job'
2020-08-22 14:05:05,065 - pyiron_log - WARNING - The job T_800 is being loaded instead of running. To re-run use the argument 'delete_existing_job=True in create_job'

The job al_temp was saved and received the ID: 51

In [47]:
al_container.structure_lst[0].info

Out[47]:
{'T': 400, 'P': 0}
In [49]:
al_container.structure_lst

Out[49]:
InputList([Al: [0.13389146 3.96541338 4.05893092]
Al: [3.99018226 2.0071096  1.95618182]
Al: [1.98560236 3.88778884 2.0465924 ]
Al: [2.04906472 2.05913422 0.09311447]
pbc: [ True  True  True]
cell:
Cell([[4.079370396328773, 2.497893949200251e-16, 2.497893949200251e-16], [0.0, 3.973148678151775, 2.4328519056175543e-16], [0.0, 0.0, 4.077409804014061]])
, Al: [0.0070279  4.03832899 0.08383998]
Al: [4.08339864 2.06533333 2.03444326]
Al: [2.20534808 4.07618808 1.94632881]
Al: [1.91118709 2.15964157 0.05514228]
pbc: [ True  True  True]
cell:
Cell([[4.103480856873615, 2.5126573483663555e-16, 2.5126573483663555e-16], [0.0, 4.113163987813141, 2.518586556021763e-16], [0.0, 0.0, 4.11975432838738]])
, Al: [3.7382874  0.12171228 4.27645154]
Al: [0.05199557 1.91099383 2.20493355]
Al: [1.92074788 0.03592662 2.13915097]
Al: [1.89264518 1.93451826 0.04368514]
pbc: [ True  True  True]
cell:
Cell([[3.801838019513017, 2.3279543807366644e-16, 2.3279543807366644e-16], [0.0, 4.003150985990484, 2.4512230207484086e-16], [0.0, 0.0, 4.332110602330073]])
])
In [ ]: