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
import numpy as np
%matplotlib inline
import matplotlib.pylab as plt
and create a pyiron project named 'structures':
from pyiron import Project
pr = Project(path='structures')
In this section we discuss various possibilities to create bulk crystal structures.
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.
structure = pr.create_structure('Al',
bravais_basis='fcc',
lattice_constant=4.05)
To plot the structure interactively in 3D simply use:
structure.plot3d()
NGLWidget()
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
structure = pr.create_ase_bulk('Al', cubic=True)
structure.plot3d()
NGLWidget()
Example: wurtzite GaN in a 3x3x3 repeated orthorhombic cell.
Note:
spacefill=False
in the plot3d()
method changes the atomic structure style to "ball and stick".structure = pr.create_ase_bulk('AlN',
crystalstructure='wurtzite',
a=3.5, orthorhombic=True)
structure.set_repeat([3,3,3])
structure.plot3d(spacefill=False)
NGLWidget()
Example: Creating a 3x4 fcc Al(111) surface with 4 layers and a vacuum of 10 Ångström
Al_111 = pr.create_surface("Al", surface_type="fcc111",
size=(3, 4, 4), vacuum=10, orthogonal=True)
Al_111.plot3d()
NGLWidget()
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.
from pyiron import create_ase_bulk, create_surface
structure = create_ase_bulk('AlN',
crystalstructure='wurtzite',
a=3.5, orthorhombic=True)
structure.set_repeat([3,3,3])
structure.plot3d(spacefill=False)
NGLWidget()
Al_111 = create_surface("Al", surface_type="fcc111",
size=(3, 4, 4), vacuum=10, orthogonal=True)
Al_111.plot3d()
NGLWidget()
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)
skutterudite.plot3d()
NGLWidget()
Using the bulk aluminum fcc example from before the structure object can be created by
structure = pr.create_ase_bulk('Al', cubic=True)
A summary of the information about the structure is given by using
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: Cell([4.05, 4.05, 4.05])
The cell vectors of the structure object can be accessed and edited through
structure.cell
Cell([4.05, 4.05, 4.05])
The positions of the atoms in the structure object can be accessed and edited through
structure.positions
array([[0. , 0. , 0. ], [0. , 2.025, 2.025], [2.025, 0. , 2.025], [2.025, 2.025, 0. ]])
We start by setting up a 4x4x4 supercell
structure = pr.create_ase_bulk('Al', cubic=True)
structure.set_repeat([4,4,4])
To create the vacancy at position index "0" simply use:
del structure[0]
To plot the structure that now contains a vacancy run:
structure.plot3d()
NGLWidget()
# 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
# 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]
# 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
NGLWidget()
# 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
n_sub = 24
structure[np.random.permutation(len(structure))[:n_sub]] = 'Mg'
# 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
NGLWidget()
You can also set-up structures through the explicit input of the cell parameters and positions
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()
NGLWidget()
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)
# 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)
cif_structure = """\
#------------------------------------------------------------------------------
#$Date: 2015-01-27 21:58:39 +0200 (Tue, 27 Jan 2015) $
#$Revision: 130149 $
#$URL: svn://www.crystallography.net/cod/cif/1/00/87/1008753.cif $
#------------------------------------------------------------------------------
#
# This file is available in the Crystallography Open Database (COD),
# http://www.crystallography.net/
#
# All data on this site have been placed in the public domain by the
# contributors.
#
data_1008753
loop_
_publ_author_name
'Buerger, M J'
'Klein, G E'
'Donnay, G'
_publ_section_title
;
Determination of the crystal structure of nepheline
;
_journal_coden_ASTM AMMIAY
_journal_name_full 'American Mineralogist'
_journal_page_first 805
_journal_page_last 818
_journal_volume 39
_journal_year 1954
_chemical_formula_structural 'K Na3 Al4 Si4 O16'
_chemical_formula_sum 'Al4 K Na3 O16 Si4'
_chemical_name_mineral Nepheline
_chemical_name_systematic 'Potassium trisodium tetraaluminium silicate'
_space_group_IT_number 173
_symmetry_cell_setting hexagonal
_symmetry_Int_Tables_number 173
_symmetry_space_group_name_Hall 'P 6c'
_symmetry_space_group_name_H-M 'P 63'
_cell_angle_alpha 90
_cell_angle_beta 90
_cell_angle_gamma 120
_cell_formula_units_Z 2
_cell_length_a 10.01
_cell_length_b 10.01
_cell_length_c 8.405
_cell_volume 729.4
_cod_database_code 1008753
loop_
_symmetry_equiv_pos_as_xyz
x,y,z
-y,x-y,z
y-x,-x,z
-x,-y,1/2+z
y,y-x,1/2+z
x-y,x,1/2+z
loop_
_atom_site_label
_atom_site_type_symbol
_atom_site_symmetry_multiplicity
_atom_site_Wyckoff_symbol
_atom_site_fract_x
_atom_site_fract_y
_atom_site_fract_z
_atom_site_occupancy
_atom_site_attached_hydrogens
_atom_site_calc_flag
K1 K1+ 2 a 0. 0. 0. 1. 0 d
Al1 Al3+ 2 b 0.3333 0.6667 0.18 1. 0 d
Si1 Si4+ 2 b 0.3333 0.6667 0.82 1. 0 d
O1 O2- 2 b 0.3333 0.6667 0. 1. 0 d
Na1 Na1+ 6 c 0.008 0.432 0. 1. 0 d
Al2 Al3+ 6 c 0.092 0.33 0.67 1. 0 d
Si2 Si4+ 6 c 0.092 0.33 0.33 1. 0 d
O2 O2- 6 c 0.02 0.33 0.5 1. 0 d
O3 O2- 6 c 0.18 0.5 0.75 1. 0 d
O4 O2- 6 c 0.17 0.53 0.25 1. 0 d
O5 O2- 6 c 0.23 0.28 0.25 1. 0 d
O6 O2- 6 c 0.23 0.28 0.75 1. 0 d
loop_
_atom_type_symbol
_atom_type_oxidation_number
K1+ 1.000
Al3+ 3.000
Si4+ 4.000
O2- -2.000
Na1+ 1.000"""
# Download and save the structure file locally
# import urllib
# urllib.request.urlretrieve(url=url, filename='strucs.'+filename);
with open('strucs.'+filename, "w") as f:
f.writelines(cif_structure)
# Using ase parsers to read the structure and then convert to a pyiron instance
import ase
from pyiron import ase_to_pyiron
structure = ase_to_pyiron(ase.io.read(filename='strucs.'+filename,
format='cif'))
structure.info["cod"] = cod
/srv/conda/envs/notebook/lib/python3.7/site-packages/ase/io/cif.py:380: UserWarning: crystal system 'hexagonal' is not interpreted for space group Spacegroup(173, setting=1). This may result in wrong setting! setting_name, spacegroup))
structure.plot3d()
NGLWidget()
Structures can be stored indepently from jobs in HDF5 by using the special StructureContainer
job. To save to disk, call run()
.
container = pr.create_job(pr.job_type.StructureContainer, "nepheline")
container.structure = structure
container.run()
The job nepheline was saved and received the ID: 1
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.
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 The job T_400 was saved and received the ID: 2 The job T_600 was saved and received the ID: 3 The job T_800 was saved and received the ID: 4 The job al_temp was saved and received the ID: 5
al_container.structure_lst[0].info
{'T': 400, 'P': 0}
al_container.structure_lst
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.077409804014059]]) , 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.103480856873612, 2.5126573483663535e-16, 2.5126573483663535e-16], [0.0, 4.11316398781314, 2.5185865560217624e-16], [0.0, 0.0, 4.119754328387385]]) , 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.8018380195130206, 2.3279543807366664e-16, 2.3279543807366664e-16], [0.0, 4.003150985990483, 2.451223020748408e-16], [0.0, 0.0, 4.332110602330072]]) ])