Introduction

This notebook was written by Janine George (E-mail: janine.george@bam.de BAM, https://jageo.github.io/).

This notebook shows how to plot Crystal Orbital Hamilton Population (COHP) and projected densities of states calculated with the Local-Orbital Basis Suite Towards Electronic-Structure Reconstruction (LOBSTER) code. Furtheremore, the classes Icohplist and Charge to evaluate ICOHPLIST.lobster and CHARGE.lobster are explained. See http://www.cohp.de for more information. The code to plot COHP and evaluate ICOHPLIST.lobster in pymatgen was started Marco Esters and extended by Janine George. The classes to plot DOSCAR.lobster, and to evaluate CHARGE.lobster were written by Janine George.

You can find our publication describing all modules and the implementation here: [1] J. George, G. Petretto, A. Naik, M. Esters, A. J. Jackson, R. Nelson, R. Dronskowski, G.-M. Rignanese, G. Hautier, ChemPlusChem n.d., n/a, e202200123.

# Uncomment the subsequent lines in this cell to install dependencies for Google Colab.
# !pip install pymatgen==2022.7.25

How to plot COHPCAR.lobster

get relevant classes

from pymatgen.electronic_structure.cohp import CompleteCohp
from pymatgen.electronic_structure.plotter import CohpPlotter

%matplotlib inline

get a completecohp object to simplify the plotting

COHPCAR_path = "lobster_data/GaAs/COHPCAR.lobster_new"
POSCAR_path = "lobster_data/GaAs/POSCAR_new"

completecohp = CompleteCohp.from_file(
    fmt="LOBSTER", filename=COHPCAR_path, structure_file=POSCAR_path
)

3809901

plot certain COHP

You have to search for the label of the COHP you would like to plot in ICOHPLIST.lobster

# search for the number of the COHP you would like to plot in ICOHPLIST.lobster (the numbers in COHPCAR.lobster are different!)
label = "16"
cp = CohpPlotter()
# get a nicer plot label
plotlabel = (
    str(completecohp.bonds[label]["sites"][0].species_string)
    + "-"
    + str(completecohp.bonds[label]["sites"][1].species_string)
)

cp.add_cohp(plotlabel, completecohp.get_cohp_by_label(label=label))
# check which COHP you are plotting

print(
    "This is a COHP between the following sites: "
    + str(completecohp.bonds[label]["sites"][0])
    + " and "
    + str(completecohp.bonds[label]["sites"][1])
)

x = cp.get_plot(integrated=False)
x.ylim([-10, 6])

x.show()

This is a COHP between the following sites: [0. 0. 0.] Ga and [0. 0. 0.] Ga

add several COHPs

# labels of the COHPs that will be summed!
labelist = ["16", "21"]
cp = CohpPlotter()
# get a nicer plot label
plotlabel = "two bonds"

cp.add_cohp(
    plotlabel,
    completecohp.get_summed_cohp_by_label_list(label_list=labelist, divisor=1),
)
x = cp.get_plot(integrated=False)
x.ylim([-10, 6])
x.show()

focus on certain orbitals only

# search for the number of the COHP you would like to plot in ICOHPLIST.lobster (the numbers in COHPCAR.lobster are different!)
label = "16"
cp = CohpPlotter()

# get orbital object
from pymatgen.electronic_structure.core import Orbital

# interaction between 4s and 4px
orbitals = [[4, Orbital.s], [4, Orbital.py]]
orbitals2 = [[4, Orbital.s], [4, Orbital.pz]]
# get a nicer plot label

plotlabel = (
    str(completecohp.bonds[label]["sites"][0].species_string)
    + "(4s)"
    + "-"
    + str(completecohp.bonds[label]["sites"][1].species_string)
    + "(4py)"
)
plotlabel2 = (
    str(completecohp.bonds[label]["sites"][0].species_string)
    + "(4s)"
    + "-"
    + str(completecohp.bonds[label]["sites"][1].species_string)
    + "(4pz)"
)

cp.add_cohp(
    plotlabel, completecohp.get_orbital_resolved_cohp(label=label, orbitals=orbitals)
)
cp.add_cohp(
    plotlabel2, completecohp.get_orbital_resolved_cohp(label=label, orbitals=orbitals2)
)
# check which COHP you are plotting

# with integrated=True, you can plot the integrated COHP
x = cp.get_plot(integrated=False)
x.ylim([-10, 6])

x.show()

How to evaluate ICOHPLIST.lobster

get relevant classes

from pymatgen.io.lobster import Icohplist

read in ICOHPLIST.lobster and get Icohpcollection object

icohplist = Icohplist(filename="lobster_data/GaAs/ICOHPLIST.lobster")
icohpcollection = icohplist.icohpcollection

get interesting properties from ICOHPLIST.lobster

# get icohp value by label (labelling according to ICOHPLIST.lobster)
# for spin polarized calculations you can also sum the spin channels
print("icohp value for certain bond by label")
label = "16"
print(icohpcollection.get_icohp_by_label(label))
print()
# you can get all Icohpvalue objects for certain bond lengths
print("Icohp values for certain bonds with certain bond lengths")
for key, icohp in icohpcollection.get_icohp_dict_by_bondlengths(
    minbondlength=0.0, maxbondlength=3.0
).items():
    print(key + ":" + str(icohp.icohp))
print()
# you can get all icohps for a certain site
print("ICOHP values of certain site")
for key, icohp in icohpcollection.get_icohp_dict_of_site(
    site=0, minbondlength=0.0, maxbondlength=3.0
).items():
    print(key + ":" + str(icohp.icohp))

icohp value for certain bond by label
-4.35464

Icohp values for certain bonds with certain bond lengths
16:{<Spin.up: 1>: -4.35464}
21:{<Spin.up: 1>: -4.35464}
23:{<Spin.up: 1>: -4.35464}
24:{<Spin.up: 1>: -4.35464}

ICOHP values of certain site
16:{<Spin.up: 1>: -4.35464}
21:{<Spin.up: 1>: -4.35464}
23:{<Spin.up: 1>: -4.35464}
24:{<Spin.up: 1>: -4.35464}

How to plot DOSCAR.lobster:

get relevant classes

from pymatgen.core.composition import Element
from pymatgen.electronic_structure.plotter import DosPlotter

# relevant classes
from pymatgen.io.lobster import Doscar

%matplotlib inline

read in DOSCAR.lobster and get structure object for later

# read in DOSCAR.lobster
doscar = Doscar(
    doscar="lobster_data/GaAs/DOSCAR.lobster", structure_file="lobster_data/GaAs/POSCAR"
)
complete_dos = doscar.completedos
# get structure object
structure = complete_dos.structure

plot total density of states

# plot total dos
Plotter = DosPlotter()
Plotter.add_dos("Total Dos", doscar.tdos)
Plotter.get_plot().show()

plot DOS projected on s, p, and d orbitals for certain element

# plot DOS of s,p, and d orbitals for certain element
Plotter = DosPlotter()
el = Element("Ga")
Plotter.add_dos_dict(complete_dos.get_element_spd_dos(el=el))
Plotter.get_plot().show()

plot DOS for cetain sites and orbitals

Plotter = DosPlotter()
# choose the sites you would like to plot
for isite, site in enumerate(structure[0:1]):
    # name the orbitals you would like to include
    # the other orbitals are named in a similar way. The orbitals are called: "s", "p_y", "p_z", "p_x", "d_xy", "d_yz", "d_z^2","d_xz", "d_x^2-y^2", "f_y(3x^2-y^2)", "f_xyz","f_yz^2", "f_z^3", "f_xz^2", "f_z(x^2-y^2)", "f_x(x^2-3y^2)"
    for orbital in ["4s"]:
        Plotter.add_dos(
            "Ga" + str(isite + 1) + ":" + orbital,
            complete_dos.get_site_orbital_dos(site, orbital),
        )
Plotter.get_plot().show()