How to use SMACT - 10 common examples

import smact.distorter as distort
import numpy as np
import ase.io as io

# Build the input
test_case = builder.cubic_perovskite(['Ba','Ti','O'],repetitions=[2,2,2])


print "------------------------------"
print "Original coordinates: ", test_case
print "------------------------------"


# Do the single substitution first, it is trivial as all Ba sites are equivalent we will choose the first Ba
subs_site = [0.0, 0.0, 0.0]  # Fractional coordinates of the site to substitute
single_substitution = distort.make_substitution(test_case,subs_site,"Sr")
print "Single: ", single_substitution


#Build a sub-lattice you wish to disorder [test case do the Ba sub-lattice]
sub_lattice = distort.build_sub_lattice(single_substitution,"Ba")


# Enumerate the inequivalent sites
inequivalent_sites = distort.get_inequivalent_sites(sub_lattice,single_substitution)


# Replace Ba at inequivalent sites with Sr
i = 0
for inequivalent_site in inequivalent_sites:
    print "------------------------------"
    print " Substituted coordinates"

#Build a sub-lattice you wish to disorder [test case do the Ba sub-lattice]
sub_lattice = distort.build_sub_lattice(single_substitution,"Ba")


# Enumerate the inequivalent sites
inequivalent_sites = distort.get_inequivalent_sites(sub_lattice,single_substitution)


# Replace Ba at inequivalent sites with Sr
i = 0
for inequivalent_site in inequivalent_sites:
    print "------------------------------"
    print " Substituted coordinates" 
    #print test_case,inequivalent_site
    distorted = distort.make_substitution(single_substitution,inequivalent_site,"Sr")
    io.write('POSCAR-%s'%i,distorted,format='vasp')
    i = i + 1
#    print distorted
    

print "------------------------------"
print "------------------------------"

import ase
from ase.io import *
import smact.builder as build
import smact.surface as surface

PbZrO = build.cubic_perovskite(['Pb','Zr','O'])

S100 = surface.cut100(PbZrO)
S110 = surface.cut110(PbZrO)
S111 = surface.cut111(PbZrO)

write('100.cif',S100)
write('110.cif',S110)
write('111.cif',S111)
write('Bulk.cif',PbZrO)

# Example script of using distorter, generate all possible (symmetry inequivalent) subsitiutions of Sr on Ba 
# sites; single and double substitutions.

import ase
import smact.builder as builder
import smact.distorter as distort
import numpy as np
import ase.io as io

# Build the input
test_case = builder.cubic_perovskite(['Ba','Ti','O'],repetitions=[2,2,2])


print "------------------------------"
print "Original coordinates: ", test_case
print "------------------------------"


# Do the single substitution first, it is trivial as all Ba sites are equivalent we will choose the first Ba
subs_site = [0.0, 0.0, 0.0]  # Fractional coordinates of the site to substitute
single_substitution = distort.make_substitution(test_case,subs_site,"Sr")
print "Single: ", single_substitution


#Build a sub-lattice you wish to disorder [test case do the Ba sub-lattice]
sub_lattice = distort.build_sub_lattice(single_substitution,"Ba")

test_case = builder.cubic_perovskite(['Ba','Ti','O'],repetitions=[2,2,2])


print "------------------------------"
print "Original coordinates: ", test_case
print "------------------------------"


# Do the single substitution first, it is trivial as all Ba sites are equivalent we will choose the first Ba
subs_site = [0.0, 0.0, 0.0]  # Fractional coordinates of the site to substitute
single_substitution = distort.make_substitution(test_case,subs_site,"Sr")
print "Single: ", single_substitution


#Build a sub-lattice you wish to disorder [test case do the Ba sub-lattice]
sub_lattice = distort.build_sub_lattice(single_substitution,"Ba")


# Enumerate the inequivalent sites
inequivalent_sites = distort.get_inequivalent_sites(sub_lattice,single_substitution)


# Replace Ba at inequivalent sites with Sr
i = 0
for inequivalent_site in inequivalent_sites:
    print "------------------------------"
    print " Substituted coordinates" 
    #print test_case,inequivalent_site
    distorted = distort.make_substitution(single_substitution,inequivalent_site,"Sr")
    io.write('POSCAR-%s'%i,distorted,format='vasp')
    i = i + 1
#    print distorted

print "Original coordinates: ", test_case
print "------------------------------"


# Do the single substitution first, it is trivial as all Ba sites are equivalent we will choose the first Ba
subs_site = [0.0, 0.0, 0.0]  # Fractional coordinates of the site to substitute
single_substitution = distort.make_substitution(test_case,subs_site,"Sr")
print "Single: ", single_substitution


#Build a sub-lattice you wish to disorder [test case do the Ba sub-lattice]
sub_lattice = distort.build_sub_lattice(single_substitution,"Ba")


# Enumerate the inequivalent sites
inequivalent_sites = distort.get_inequivalent_sites(sub_lattice,single_substitution)


# Replace Ba at inequivalent sites with Sr
i = 0
for inequivalent_site in inequivalent_sites:
    print "------------------------------"
    print " Substituted coordinates" 
    #print test_case,inequivalent_site
    distorted = distort.make_substitution(single_substitution,inequivalent_site,"Sr")
    io.write('POSCAR-%s'%i,distorted,format='vasp')
    i = i + 1
#    print distorted
    

print "------------------------------"
print "------------------------------"

count_of_ternary = 0
count_raw_ternary = 0
element_count = 0
ion_count = 0
charge_neutral_count = 0
pauling_count = 0
for i, ele_a in enumerate(list_of_elements):
 	element_count = element_count + 1
	paul_a = smact.Element(ele_a).pauling_eneg
	for ox_a in smact.Element(ele_a).oxidation_states:
		ion_count = ion_count + 1
		for j, ele_b in enumerate(list_of_elements[i+1:]):
			element_count = element_count + 1
			paul_b = smact.Element(ele_b).pauling_eneg
			for ox_b in smact.Element(ele_b).oxidation_states:
				ion_count = ion_count + 1
				element = [ele_a, ele_b]
				print element
				oxidation_states = [ox_a, ox_b]
				pauling_electro = [paul_a, paul_b]
				electroneg_makes_sense = pauling_test(oxidation_states, pauling_electro)
				cn_e, cn_r = smact.charge_neutrality([ox_a, ox_b],threshold=1)
				if cn_e:
					charge_neutral_count = charge_neutral_count + 1
					if electroneg_makes_sense:
						pauling_count = pauling_count + 1

print "  "
print "------------------------------------------"
print "------------------------------------------"
print "Summary of screening for ternary compounds"

species = line.split()
		if int(species[0]) > 0 and int(species[0]) < 81:
			#if len(smact.Element(species[1]).oxidation_states) > 0:
			if smact.Element(species[1]).oxidation_states:
			    list_of_elements.append(species[1])
			
count_of_ternary = 0
count_raw_ternary = 0
element_count = 0
ion_count = 0
charge_neutral_count = 0
pauling_count = 0
for i, ele_a in enumerate(list_of_elements):
 	element_count = element_count + 1
	paul_a = smact.Element(ele_a).pauling_eneg
	for ox_a in smact.Element(ele_a).oxidation_states:
		ion_count = ion_count + 1
		for j, ele_b in enumerate(list_of_elements[i+1:]):
			element_count = element_count + 1
			paul_b = smact.Element(ele_b).pauling_eneg
			for ox_b in smact.Element(ele_b).oxidation_states:
				ion_count = ion_count + 1
				for k, ele_c in enumerate(list_of_elements[i+j+2:]):
				  	element_count = element_count + 1
					paul_c = smact.Element(ele_c).pauling_eneg
					for ox_c in smact.Element(ele_c).oxidation_states:
						ion_count = ion_count + 1
						element = [ele_a, ele_b, ele_c]
						print element
						oxidation_states = [ox_a, ox_b, ox_c]
						pauling_electro = [paul_a, paul_b, paul_c]
						electroneg_makes_sense = pauling_test(oxidation_states, pauling_electro)

list_of_elements.append(species[1])
			
count_of_ternary = 0
count_raw_ternary = 0
element_count = 0
ion_count = 0
charge_neutral_count = 0
pauling_count = 0
for i, ele_a in enumerate(list_of_elements):
 	element_count = element_count + 1
	paul_a = smact.Element(ele_a).pauling_eneg
	for ox_a in smact.Element(ele_a).oxidation_states:
		ion_count = ion_count + 1
		for j, ele_b in enumerate(list_of_elements[i+1:]):
			element_count = element_count + 1
			paul_b = smact.Element(ele_b).pauling_eneg
			for ox_b in smact.Element(ele_b).oxidation_states:
				ion_count = ion_count + 1
				for k, ele_c in enumerate(list_of_elements[i+j+2:]):
				  	element_count = element_count + 1
					paul_c = smact.Element(ele_c).pauling_eneg
					for ox_c in smact.Element(ele_c).oxidation_states:
						ion_count = ion_count + 1
						element = [ele_a, ele_b, ele_c]
						print element
						oxidation_states = [ox_a, ox_b, ox_c]
						pauling_electro = [paul_a, paul_b, paul_c]
						electroneg_makes_sense = pauling_test(oxidation_states, pauling_electro)
						cn_e, cn_r = smact.charge_neutrality([ox_a, ox_b, ox_c],threshold=1)
						if cn_e:
							charge_neutral_count = charge_neutral_count + 1
							if electroneg_makes_sense:

list_of_elements.append(species[1])
			
count_of_ternary = 0
count_raw_ternary = 0
element_count = 0
ion_count = 0
charge_neutral_count = 0
pauling_count = 0
for i, ele_a in enumerate(list_of_elements):
 	element_count = element_count + 1
	paul_a = smact.Element(ele_a).pauling_eneg
	for ox_a in smact.Element(ele_a).oxidation_states:
		ion_count = ion_count + 1
		for j, ele_b in enumerate(list_of_elements[i+1:]):
			element_count = element_count + 1
			paul_b = smact.Element(ele_b).pauling_eneg
			for ox_b in smact.Element(ele_b).oxidation_states:
				ion_count = ion_count + 1
				element = [ele_a, ele_b]
				print element
				oxidation_states = [ox_a, ox_b]
				pauling_electro = [paul_a, paul_b]
				electroneg_makes_sense = pauling_test(oxidation_states, pauling_electro)
				cn_e, cn_r = smact.charge_neutrality([ox_a, ox_b],threshold=1)
				if cn_e:
					charge_neutral_count = charge_neutral_count + 1
					if electroneg_makes_sense:
						pauling_count = pauling_count + 1

print "  "
print "------------------------------------------"
print "------------------------------------------"