#!/usr/bin/env python
# coding: utf-8

# # 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](https://wiki.fysik.dtu.dk/ase/)). This makes it possible to use routines from ASE to help set-up structures.
# 
# Furthermore, pyiron uses the [NGLview](http://nglviewer.org/nglview/latest/api.html) 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
get_ipython().run_line_magic('matplotlib', 'inline')
import matplotlib.pylab as plt


# and create a pyiron project named 'structures':

# In[2]:


from pyiron 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](https://wiki.fysik.dtu.dk/ase/ase/build/build.html#ase.build.bulk). 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](https://wiki.fysik.dtu.dk/ase/_modules/ase/build/surface.html)

# **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)


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

# In[15]:


structure.cell


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

# In[16]:


structure.positions


# ## 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())


# 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()


# ### 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()


# ## 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)](http://www.crystallography.net/cod/index.php)

# 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]:


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"""


# In[29]:


# 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)


# In[30]:


# 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


# In[31]:


structure.plot3d()


# 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()


# 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[33]:


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()


# In[34]:


al_container.structure_lst[0].info


# In[35]:


al_container.structure_lst


# In[ ]: