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.

In [1]:
# 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¶

In [2]:
from pymatgen.electronic_structure.cohp import CompleteCohp
from pymatgen.electronic_structure.plotter import CohpPlotter

%matplotlib inline

get a completecohp object to simplify the plotting¶

In [3]:
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

In [4]:
# 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¶

In [5]:
# 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¶

In [6]:
# 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¶

In [7]:
from pymatgen.io.lobster import Icohplist

read in ICOHPLIST.lobster and get Icohpcollection object¶

In [8]:
icohplist = Icohplist(filename="lobster_data/GaAs/ICOHPLIST.lobster")
icohpcollection = icohplist.icohpcollection

get interesting properties from ICOHPLIST.lobster¶

In [9]:
# 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¶

In [10]:
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¶

In [11]:
# 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¶

In [12]:
# 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¶

In [13]:
# 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¶

In [14]:
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()

evaluate CHARGE.lobster¶

get relevant classes¶

In [15]:
from pymatgen.io.lobster import Charge

read in charge and produce a structure with the charge as a property¶

In [16]:
charge = Charge(filename="lobster_data/GaAs/CHARGE.lobster")
newstructure = charge.get_structure_with_charges(
    structure_filename="lobster_data/GaAs/POSCAR"
)
print(newstructure)

Full Formula (Ga1 As1)
Reduced Formula: GaAs
abc   :   4.065993   4.065993   4.065993
angles:  60.000000  60.000000  60.000000
pbc   :       True       True       True
Sites (2)
  #  SP       a     b     c    Loewdin Charges    Mulliken Charges
---  ----  ----  ----  ----  -----------------  ------------------
  0  Ga    0     0     0                 -0.06                0.15
  1  As    0.75  0.75  0.75               0.06               -0.15

evaluate GROSSPOP.lobster¶

get relevant class¶

In [17]:
from pymatgen.io.lobster import Grosspop
In [18]:
grosspop = Grosspop(filename="lobster_data/GaAs/GROSSPOP.lobster")
print(grosspop.list_dict_grosspop)

[{'element': 'Ga', 'Mulliken GP': {'4s': 1.22, '4p_y': 0.54, '4p_z': 0.54, '4p_x': 0.54, '3d_xy': 2.0, '3d_yz': 2.0, '3d_z^2': 2.0, '3d_xz': 2.0, '3d_x^2-y^2': 2.0, 'total': 12.85}, 'Loewdin GP': {'4s': 1.1, '4p_y': 0.66, '4p_z': 0.66, '4p_x': 0.66, '3d_xy': 2.0, '3d_yz': 2.0, '3d_z^2': 2.0, '3d_xz': 2.0, '3d_x^2-y^2': 2.0, 'total': 13.06}}, {'element': 'As', 'Mulliken GP': {'4s': 1.64, '4p_y': 1.17, '4p_z': 1.17, '4p_x': 1.17, 'total': 5.15}, 'Loewdin GP': {'4s': 1.38, '4p_y': 1.19, '4p_z': 1.19, '4p_x': 1.19, 'total': 4.94}}]

get a structure with total gross populations¶

In [19]:
newstructure = grosspop.get_structure_with_total_grosspop("lobster_data/GaAs/POSCAR")
print("Structure:")
print(newstructure)

Structure:
Full Formula (Ga1 As1)
Reduced Formula: GaAs
abc   :   4.065993   4.065993   4.065993
angles:  60.000000  60.000000  60.000000
pbc   :       True       True       True
Sites (2)
  #  SP       a     b     c    Total Loewdin GP    Total Mulliken GP
---  ----  ----  ----  ----  ------------------  -------------------
  0  Ga    0     0     0                  13.06                12.85
  1  As    0.75  0.75  0.75                4.94                 5.15

FATBAND plot¶

get relevant classes¶

In [20]:
from pymatgen.electronic_structure.plotter import (
    BSDOSPlotter,
    BSPlotter,
    BSPlotterProjected,
)
from pymatgen.io.lobster import Fatband

get a bandstructure plot that is combined with a DOS plot¶

In [21]:
fatband = Fatband(
    filenames="lobster_data/GaAs",
    vasprun="lobster_data/GaAs/vasprun.xml",
    Kpointsfile="lobster_data/GaAs/KPOINTS",
)
# get a band structure object
bssymline = fatband.get_bandstructure()
# print(bssymline.as_dict())

# this can be plotted with the classes to plot bandstructures from vasp
BSDOSPlotter(bs_projection="elements", dos_projection="elements").get_plot(
    bs=bssymline, dos=complete_dos
).show()

/hpc-user/jgeorge/PycharmProjects/UpdateTutorials/pymatgen/pymatgen/io/lobster/outputs.py:1089: UserWarning: Make sure all relevant FATBAND files were generated and read in!
  warnings.warn("Make sure all relevant FATBAND files were generated and read in!")
/hpc-user/jgeorge/PycharmProjects/UpdateTutorials/pymatgen/pymatgen/io/lobster/outputs.py:1090: UserWarning: Use Lobster 3.2.0 or newer for fatband calculations!
  warnings.warn("Use Lobster 3.2.0 or newer for fatband calculations!")
findfont: Font family ['Times New Roman'] not found. Falling back to DejaVu Sans.
findfont: Font family ['Times New Roman'] not found. Falling back to DejaVu Sans.
findfont: Font family ['Times New Roman'] not found. Falling back to DejaVu Sans.
In [22]:
# another plot type from pymatgen:
bsplotter = BSPlotterProjected(bssymline)
bsplotter.get_projected_plots_dots(
    {"Ga": ["4s", "4p", "3d"], "As": ["4s", "4p"]}
).show()

Read lobsterout¶

In [23]:
from pymatgen.io.lobster import Lobsterout

get all relevant infos from lobsterout¶

In [24]:
lobsterout = Lobsterout("lobster_data/GaAs/lobsterout")

document = lobsterout.get_doc()

charge spilling can be accessed easily¶

In [25]:
print(document["chargespilling"])

[0.009899999999999999]

Create input files for vasp and lobster automatically¶

relevant classes¶

In [26]:
from pymatgen.io.lobster import Lobsterin

a Lobsterin object with standard settings is created, a standard basis is used¶

In [27]:
lobsterin = Lobsterin.standard_calculations_from_vasp_files(
    "lobster_data/GaAs/POSCAR",
    "lobster_data/GaAs/INCAR",
    "lobster_data/GaAs/POTCAR",
    option="standard",
)

/hpc-user/jgeorge/PycharmProjects/UpdateTutorials/pymatgen/pymatgen/io/lobster/inputs.py:670: UserWarning: Always check and test the provided basis functions. The spilling of your Lobster calculation might help
  warnings.warn(
/hpc-user/jgeorge/PycharmProjects/UpdateTutorials/pymatgen/pymatgen/io/vasp/inputs.py:1822: UnknownPotcarWarning: POTCAR with symbol Ga_d has metadata that does not match                          any VASP POTCAR known to pymatgen. The data in this                          POTCAR is known to match the following functionals:                          ['PBE_52', 'PBE_54', 'unvie_PBE_52']
  warnings.warn(
/hpc-user/jgeorge/PycharmProjects/UpdateTutorials/pymatgen/pymatgen/io/vasp/inputs.py:1822: UnknownPotcarWarning: POTCAR with symbol As has metadata that does not match                          any VASP POTCAR known to pymatgen. The data in this                          POTCAR is known to match the following functionals:                          ['PBE_52', 'PBE_54', 'unvie_PBE_52']
  warnings.warn(

writes lobsterin¶

In [28]:
lobsterin.write_lobsterin(path="lobsterin")
file = open("./lobsterin")
print(file.read())

COHPstartEnergy -35.0
COHPendEnergy 5.0
gaussianSmearingWidth 0.1
basisSet pbeVaspFit2015
cohpGenerator from 0.1 to 6.0 orbitalwise
saveProjectionToFile
basisfunctions Ga 3d 4p 4s 
basisfunctions As 4p 4s

will change ISYM to -1, NSW to 0, insert NBANDS, and set LWAVE to True in the INCAR¶

additional changes to the INCAR are allowed via further_settings¶

In [29]:
lobsterin.write_INCAR(
    incar_input="lobster_data/GaAs/INCAR",
    incar_output="INCAR.lobster",
    poscar_input="lobster_data/GaAs/POSCAR",
    isym=-1,
    further_settings={"IBRION": -1},
)
file = open("./INCAR.lobster")
print(file.read())

ALGO = Normal
EDIFF = 1e-07
EDIFFG = 1e-05
ENCUT = 520
IBRION = -1
ISIF = 3
ISMEAR = 0
ISPIN = 1
ISYM = -1
LCHARG = False
LORBIT = 11
LWAVE = True
NBANDS = 13
NELM = 200
NELMIN = 4
NPAR = 4
NSW = 0
POTIM = 0.5
PREC = Accurate
SIGMA = 0.1


/hpc-user/jgeorge/PycharmProjects/UpdateTutorials/pymatgen/pymatgen/io/lobster/inputs.py:314: UserWarning: Please check your incar_input before using it. This method only changes three settings!
  warnings.warn("Please check your incar_input before using it. This method only changes three settings!")

a Lobsterin object with standard settings is created, a basis given by the user is used¶

In [30]:
lobsterin = Lobsterin.standard_calculations_from_vasp_files(
    "lobster_data/GaAs/POSCAR",
    "lobster_data/GaAs/INCAR",
    "lobster_data/GaAs/POTCAR",
    option="standard",
    dict_for_basis={"Ga": "4s 4p", "As": "4s 4p"},
)

# writes lobsterin
lobsterin.write_lobsterin(path="lobsterin")
file = open("./lobsterin")
print(file.read())

COHPstartEnergy -35.0
COHPendEnergy 5.0
gaussianSmearingWidth 0.1
basisSet pbeVaspFit2015
cohpGenerator from 0.1 to 6.0 orbitalwise
saveProjectionToFile
basisfunctions Ga 4s 4p
basisfunctions As 4s 4p