In [1]:

    

from pymatgen import Structure, Molecule
from monty.json import MontyDecoder
from monty.serialization import loadfn


    

In [2]:

    

#define the coulomb energy model
from pymatgen.analysis.energy_models import EnergyModel

class ColoumbEnergy(EnergyModel):
    def get_energy(self, structure):
        """
        use the atomic numbers of the structures's elements to compute 
        the coloumbic energy of the structure
        """
        energy = 0
        for i, sitei in enumerate(structure):
            for j, sitej in enumerate(structure):
                if i != j:
                    dij = structure.get_distance(i,j)
                    Zi = sitei.species_and_occu.items()[0][0].Z
                    Zj = sitej.species_and_occu.items()[0][0].Z
                    energy += 0.5 * Zi*Zj/dij
        return energy
    
    def as_dict(self):
        return {"version": __version__,
                "@module": self.__class__.__module__,
                "@class": self.__class__.__name__,
                "init_args": {"name":"coloumb"} }


    

In [3]:

    

#examples:
acetic_acid = Molecule.from_file('../test/acetic_acid.xyz')
cscl = loadfn('/home/matk/Software/pymatgen/pymatgen/util/structures/CsCl.json', cls=MontyDecoder)
ce = ColoumbEnergy()
print ce.get_energy(acetic_acid)
print ce.get_energy(cscl)


    

    




227.515696545
256.508672698

In [4]:

    

#sort structures based on their coloumbic energies
all_structures = dict([(acetic_acid.composition.formula, ce.get_energy(acetic_acid)),
                       (cscl.composition.formula, ce.get_energy(cscl))
                       ])
sorted(all_structures.items(), key=lambda d: d[1])


    

    
Out[4]:





[(u'H4 C2 O2', 227.51569654479172), (u'Cs1 Cl1', 256.5086726976955)]

In [5]:

    

#add oxidation states to molecule

from pymatgen.core.sites import Site
from pymatgen.core.periodic_table import Specie

def add_oxidation_states(oxidation_states, molecule):
    new_sites = []
    for i, site in enumerate(molecule):
        new_sp = {}
        for el, occu in site.species_and_occu.items():
            sym = el.symbol
            new_sp[Specie(sym, oxidation_states[sym])] = occu
            new_site = Site(new_sp, site.coords,properties=site.properties)
            new_sites.append(new_site)
    return Molecule.from_sites(new_sites)


    

In [6]:

    

oxidation_states = {'C':4, 'O':-2, 'H':1}
acetic_acid_decorated = add_oxidation_states(oxidation_states, acetic_acid)
print acetic_acid_decorated


    

    




Molecule Summary (H4 C2 O2)
Reduced Formula: H2CO
Charge = 0, Spin Mult = 1
Sites (8)
1 C4+     1.079784    -0.892321     0.000000
2 C4+     0.000000     0.155304     0.000000
3 O2-     0.167514     1.358238     0.000000
4 H+     2.060132    -0.408887     0.000000
5 H+     0.977098    -1.536687     0.883544
6 H+     0.977098    -1.536687    -0.883544
7 O2-    -1.244586    -0.412073     0.000000
8 H+    -1.876454     0.335039     0.000000

In [7]:

    

#add oxidation states to structures
oxidation_states = {'Cs':1, 'Cl':-1}

cscl.add_oxidation_state_by_element(oxidation_states)

#access oxidation state and ionic radius
for site in cscl:
    specie = site.species_and_occu.keys()[0] 
    print specie.element, specie.oxi_state, specie.ionic_radius
    
#    for sp, amt in site.species_and_occu.items():
#        print sp.Z, sp.X, sp.oxidation_states, sp.ionic_radii, sp.common_oxidation_states, sp.full_electronic_structure


    

    




Cs 1 1.81 ang
Cl -1 1.67 ang

In [8]:

    

###oxidation state decorations using pymatgen.core.transformations
##manual oxiadtion state decoration
#only works for Structure objects
#from pymatgen.transformations.standard_transformations import OxidationStateDecorationTransformation
#oxdec = OxidationStateDecorationTransformation({'Cs':1, 'Cl':-1})
#s_ox = oxdec.apply_transformation(cscl)

##automatic oxidation state decoration
##only works for Structure objects
#from pymatgen.transformations.standard_transformations import AutoOxiStateDecorationTransformation
#auto_oxdec = AutoOxiStateDecorationTransformation()
#print auto_oxdec.apply_transformation(cscl)


    

In [8]: