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

# ## **crystal_torture:** A Python/Fortran crystal tortuosity module  
#  `crystal_torture` is a Python, Fortran and OpenMP crystal structure analysis module. The module contains an interface to [pymatgen](http://www.pymatgen.org) and a set of classes that enables:
# 
# * a crystal structure to be converted into a simple graph for network analysis
# * connected clusters of crystal sites (nodes) to be retrieved and output
# * periodicity of connected clusters of crystal sites to be determined
# * relative path tortuosity to traverse a crystal within a connected cluster to be calculated for each site in the cluster
# 
# The module is written in Python and Python wrapped Fortran (with OpenMP).
#                                                                        

# ### **Crystals as Graphs, Connectivity and Relative Tortuosity**
# 
# A crystal structure can be described as a graph, consisting of clusters, which in turn are comprised of connected nodes. Taking a crystal structure we can treat each site within it it as a node. 
# 
# <img src="Images/crystal.png" alt="drawing" width="200px" align="middle"/>
# 
# Nodes can be considered connected or disconnected based on some user determined inter-nodal distance and occupying species criteria.
# 
# Selecting only <img style="float: margin: 0px 0px 0px 0px" ;src="Images/node.png" width="30px"/> nodes:
# 
# <img src="Images/crystal_node.png" alt="drawing" width="200px" align="middle"/>
# 
# we obtain a graph comprised of connected clusters that extends across the periodic crystal boundary:
# 
# <img src="Images/crystal_clusters.png" alt="drawing" width="200px" align="middle"/>
# 
# Selecting only the periodic clusters:
# 
# <img src="Images/crystal_cluster_periodic.png" alt="drawing" width="251px"  margin-left=-100px align="middle"/>
# 
# We can select one node and calculate the minimum number of internodal steps to traverse the full crystal and return to a periodic nodal site:
# 
# <img src="Images/crystal_cluster_periodic_pathway.png" alt="drawing" width="251px" align="middle"/>
# 
# This is what we define as the site tortuosity, &tau;. The relative site tortuosity therefore is the ratio of this value to the minimum pathway if all nodes were available to form the pathway. 
# 

# ### **Setting up a graph and torturing it**
# 
# Using pymatgen we can set up a graph structure directly from a <a href="http://pymatgen.org">pymatgen</a> compatible <a href="http://pymatgen.org/_modules/pymatgen/core/structure.html#IMolecule.from_file">filetype</a>  (in this instance a <a href="http://cms.mpi.univie.ac.at/vasp/guide/node59.html">POSCAR</a> file). We have stipulated an internodal distance of 4.0 &#8491; and selected only sites occupied by Li atoms as nodes:

# In[1]:


import sys
import time
sys.path.append('../crystal_torture')

from crystal_torture.pymatgen_interface import graph_from_file

graph = graph_from_file("POSCAR.vasp",4.0,["Li"])


# In[2]:


print(len(graph.clusters))


# We can list the indices of nodes within each cluster:

# In[3]:


for cluster in graph.clusters:
    print("Cluster: (len) (indices in cluster))")
    print( cluster.return_uc_indices())


# We can also list the periodicity of the clusters within the graph (i.e. whether the cluster traverses the entire unit cell in 1,2 or 3 dimensions):

# In[4]:


print([cluster.periodic for cluster in graph.clusters])


# We can print the fraction of nodes in the graph that are in a percolating (i.e. periodic) cluster:

# In[5]:


print('{} {}\n'.format("Frac Perc:",graph.return_frac_percolating()))


# We can torture the periodic clusters in the graph, meaning for each node in a cluster we can calculate the tortuosity of the minimum pathway from the node to any of its periodic images (i.e. the number of internodal steps to reach a periodic image). The ratio of this tortuosity to the minimum pathway is the relative site tortuosity. 
# 
# The algorithm performed is a <a href="https://en.wikipedia.org/wiki/Breadth-first_search"> breadth first search</a>. Depending on system size, this algorithm can be computationally expensive and there are two ways of running it. There is a pure python version, which can be used for small systems:
# 

# In[6]:


start_time = time.time()
graph.torture_py()
end_time = time.time()

print("Time:", end_time - start_time)        


# Or a version that performs the search using Fortran90, coupled with OpenMP. The algorithm in this case is 
# parallelised across nodes in the graph, calculating the tortuosity for each node in parallel, and therefore will perform the search much more quickly for large systems.

# In[7]:


graph = graph_from_file("POSCAR.vasp",4.0,["Li"])

time1=time.time()

start_time = time.time()
graph.torture()
end_time = time.time()

print("Time:", end_time - start_time)
        


# Following the tortuosity analysis we can see the average tortuosity for the nodes in a cluster:

# In[8]:


for no,cluster in enumerate(graph.minimal_clusters):
        print('{}{}{}\n'.format("****** Cluster: ",no, " ********"))
        print('{} {}\n'.format("Size:",cluster.size))
        print('{} {}\n'.format("Tortuosity:",cluster.tortuosity))
        print('{} {}\n'.format("Periodicity:",cluster.periodic))



# Or indeed exmaine the tortuosity of every node in a cluster:

# In[9]:


print(graph.tortuosity)
    
    


# It is also possible to output the clusters found to a filetype which can be viewed in a crystal structure viewer such as <a href="http://www.ks.uiuc.edu/Research/vmd/">VMD</a> or <a href="http://jp-minerals.org/vesta/en/">VESTA</a>

# In[10]:


graph.output_clusters(fmt='poscar',periodic=False)


# ### **Worked example: Doping and torturing Spinel MgAl<sub>2</sub>O<sub>4</sub>**
# 
# Taking the spinel structure, there are 56 ions in the conventional unit cell with stoichiometry (A<sup>2+</sup>)[B<sup>3+</sup>]O<sup>4+</sup>. In normal magnesium spinel (Mg)[Al<sub>2</sub>]O<sub>4</sub>, the Mg<sup>2+</sup> reside on the (A) sites, Al<sup>3+</sup> reside on the [B]-sites.
# 
# Taking this spinel structured magnesium spinel we can co-dope with {Li-Al}, upon which pairs of Mg<sup>2+</sup> cations are substituted in equal proportions by Li<sup>+</sup> and Al<sup>3+</sup>, giving
#     a composition of (Al<sub>x</sub>Mg<sub>1-2x</sub>Li<sub>x</sub>)[Al<sub>2</sub>]O<sub>4</sub>
# 
# <div align="center">
#     <table cellpadding="100" cellspacing="200" background-colour="white">
#         <tr>
#             <td> <div style="font-size:120%; text-align:center;"><img src="Images/POSCAR_undoped.png" alt="alternate text" height=200  style="padding-bottom:0.5cm;"/> MgAl<sub>2</sub>O<sub>4</sub></div> </td>
#             <td> <div style="font-size:120%; text-align:center;"><img src="Images/Rmoved.png" alt="alternate text" height=200 style="padding-bottom:0.5cm;"/>Doping removes two Mg<sup>s+</sup> and inserts one Li<sup>+</sup> and one Al<sup>3+</sup></div> </td>
#             <td> <div style="font-size:120%; text-align:center;"><img src="Images/Doped.png" alt="alternate text" height=200 style="padding-bottom:0.5cm;"/>Doped structure (x=0.25)(Li<sub>x</sub>Mg<sub>1−2x</sub>Al<sub>2+x</sub>O<sub>4</sub>)</div> </td>
#         </tr>
#     </table>
# </div>    
# 
# One potential use for this material could be as a solid Li<sup>+</sup> electrolyte in Li-ion batteries. 
# 
# In this case we are interested in the connectivity of the Li-ions within the crystal, and the tortuosity of possible pathways through the crystal for these Li-ions. We wish to examine the networks formed by Li-ions occupying the [A] sites.

# Taking the spinel unit cell we can make a pymatgen structure object:

# In[11]:


from pymatgen import Structure
spinel=Structure.from_file("POSCAR_SPINEL.vasp")


# We can label the sites of the structure:

# In[12]:


spinel.add_site_property("label",["A"]*8+["B"]*16+["O"]*32)


# and make a (3x3x3) supercell:

# In[13]:


spinel.make_supercell([3,3,3])


# and we can randomly dope the structure using crystal torture (we have included a very simple module for this task):

# In[14]:


import crystal_torture.pymatgen_doping as pd

spinel = pd.dope_structure(spinel,conc=0.9,species_to_rem="Mg",species_to_insert=["Li","Al"],label_to_remove="A")

print(spinel)


# we can output the full structure for inspection:

# In[15]:


print(spinel)
spinel.to(filename="POSCAR_full.vasp")


# now we can get a graph from the structure, including only [A] sites occuptied by Li ions:

# In[16]:


from crystal_torture.pymatgen_interface import graph_from_structure
graph = graph_from_structure(structure=spinel,rcut=4.0,elements={"Li"})


# We can torture the graph:

# In[17]:


graph.torture()


# output the clusters:

# In[18]:


graph.output_clusters('poscar',periodic=False)


# Now we can examine the size, periodicity and tortuosity of the clusters in graph:

# In[19]:


for cluster in graph.minimal_clusters:
    print("**************")
    print("Cluster size",cluster.size)
    print("Cluster Tortuosity",cluster.tortuosity)
    print("Cluster Periodicity",cluster.periodic)

print(graph.return_frac_percolating())


# In[ ]: