# useful to autoreload the module without restarting the kernel
%load_ext autoreload
%autoreload 2

from mppi import Parsers as P
import numpy as np

Tutorial of the PwParser class

This tutorial describes the usage of the class PwParser of mppi used to extract information for the XML output file data-file-schema produced by pw.

Parse of a output file of Silicon

Follow the tutorial for the QeCalculator to produce the xml file used by the PwParser.

The class is initialized by specifying the name of the xml file including its relative path

results = P.PwParser('QeCalculator_test/ecut_50.save/data-file-schema.xml')

Parse file : QeCalculator_test/ecut_50.save/data-file-schema.xml

When the object is initialized several attributes are set, for instance

results.units

'Hartree atomic units'

print(results.natoms)
print(results.atomic_species)
print(results.atomic_positions)
print(results.nbands)
print(results.nkpoints)
print(results.spin_degen)

2
{'Si': ['2.808600000000000e1', 'Si.pbe-mt_fhi.UPF']}
[['Si', [-1.2875, 1.2875, 1.2875]], ['Si', [1.2875, -1.2875, -1.2875]]]
4
32
2

We can extract some lattice-related properties. For instance

print(results.alat) # the lattice parameter in a.u.
lat = results.lattice # the array with the lattice vectors
a1,a2,a3 = lat
print(a1,a2,a3)
print(results.eval_lattice_volume()) # lattice volume in a.u.

10.3
[-5.15  0.    5.15] [0.   5.15 5.15] [-5.15  5.15  0.  ]
273.1817500000001

We can check that, since the lattice is fcc the volume is \(a_{lat}^3/4\)

results.alat**3/4

273.1817500000001

If we need the lattice vectors in units of the lattice constant the get_lattice method can be used

results.get_lattice(rescale=True)

array([[-0.5,  0. ,  0.5],
       [ 0. ,  0.5,  0.5],
       [-0.5,  0.5,  0. ]])

We can also compute the vectors of the reciprocal lattice

rec_lat = results.get_reciprocal_lattice(rescale=True)
b1,b2,b3 = rec_lat
print(b1,b2,b3)

[-1. -1.  1.] [1. 1. 1.] [-1.  1. -1.]

We use the tools of the ParsersUtils module to compute the lattice volume

P.ParsersUtils.eval_lattice_volume(rec_lat)

3.9999999999999987

(2*np.pi)**3/results.eval_lattice_volume()

0.9080043357303277

We check that if rescale=False, the volume of the reciprocal lattice is \((2\pi)^3/vol_{lat}\). Instead, if the option rescale = True is chosen the volume of the reciprocal lattice is 4, as expected.

We can access to the lattice symmetries

syms = results.syms
print('number of symmetries',len(syms))
syms[13]

number of symmetries 48

array([[ 0.,  0.,  1.],
       [-1., -1., -1.],
       [ 1.,  0.,  0.]])

We can print the ks energies and weight for each kpoint

for k,e,w in zip(results.kpoints,results.evals,results.weights):
    print(k,e,w)

[0. 0. 0.] [-0.21239769  0.2248015   0.2248019   0.2248019 ] [0.00925926]
[-0.16666667  0.16666667 -0.16666667] [-0.19923027  0.13749243  0.20843608  0.2084364 ] [0.05555556]
[-0.33333333  0.33333333 -0.33333333] [-0.16224325  0.02984914  0.18835302  0.18835333] [0.05555556]
[ 0.5 -0.5  0.5] [-0.12761605 -0.0295912   0.1810951   0.18109542] [0.02777778]
[0.         0.33333333 0.        ] [-0.19474764  0.14929523  0.18151447  0.18151464] [0.05555556]
[-0.16666667  0.5        -0.16666667] [-0.16503928  0.06221564  0.15212893  0.16133873] [0.11111111]
[ 0.66666667 -0.33333333  0.66666667] [-0.12174836 -0.01658335  0.12355478  0.16017609] [0.11111111]
[ 0.5        -0.16666667  0.5       ] [-0.1357136   0.00516689  0.11382729  0.17756576] [0.11111111]
[3.33333333e-01 2.77555756e-17 3.33333333e-01] [-0.17779262  0.08837439  0.13721986  0.20417699] [0.05555556]
[0.         0.66666667 0.        ] [-0.14298339  0.04081044  0.13688804  0.13688817] [0.05555556]
[ 0.83333333 -0.16666667  0.83333333] [-0.10062076 -0.01579602  0.0928995   0.13189133] [0.11111111]
[ 6.66666667e-01 -5.55111512e-17  6.66666667e-01] [-0.09359008 -0.02315423  0.06458065  0.14797813] [0.05555556]
[ 0. -1.  0.] [-0.06137332 -0.06137309  0.12102867  0.12102879] [0.02777778]
[ 0.66666667 -0.33333333  1.        ] [-0.12899013  0.01785089  0.09513332  0.14145473] [0.11111111]
[ 0.5        -0.16666667  0.83333333] [-0.09323708 -0.02570101  0.07904937  0.12455705] [0.11111111]
[-0.33333333 -1.          0.        ] [-0.05709944 -0.05709919  0.09247379  0.09247395] [0.11111111]
[-0.16666667  0.16666667  0.16666667] [-0.19923026  0.13749218  0.20843636  0.20843636] [0.01851852]
[-0.33333333  0.33333333  0.33333333] [-0.16224319  0.02984888  0.18835326  0.18835326] [0.01851852]
[ 0.5 -0.5 -0.5] [-0.12761587 -0.02959158  0.18109533  0.18109533] [0.00925926]
[ 0.16666667 -0.16666667  0.5       ] [-0.16503929  0.06221573  0.15212865  0.16133889] [0.05555556]
[-0.16666667  0.16666667  0.5       ] [-0.16503926  0.06221543  0.15212897  0.16133886] [0.05555556]
[-0.66666667  0.66666667 -0.33333333] [-0.12174836 -0.01658332  0.12355456  0.16017627] [0.05555556]
[ 0.66666667 -0.66666667 -0.33333333] [-0.12174821 -0.01658369  0.12355481  0.16017623] [0.05555556]
[-0.5         0.5        -0.16666667] [-0.13571357  0.00516682  0.11382715  0.17756596] [0.05555556]
[ 0.5        -0.5        -0.16666667] [-0.1357135   0.00516657  0.11382734  0.17756593] [0.05555556]
[-0.33333333  0.33333333  0.        ] [-0.1777926   0.08837415  0.13721989  0.20417721] [0.05555556]
[-0.83333333  0.83333333 -0.16666667] [-0.10062077 -0.01579597  0.09289929  0.13189146] [0.05555556]
[ 0.83333333 -0.83333333 -0.16666667] [-0.10062068 -0.01579618  0.09289944  0.13189145] [0.05555556]
[-0.66666667  0.66666667  0.        ] [-0.09358994 -0.02315446  0.06458057  0.14797828] [0.05555556]
[-0.66666667  1.         -0.33333333] [-0.12899018  0.0178511   0.09513332  0.14145455] [0.11111111]
[0.5        0.83333333 0.16666667] [-0.09323712 -0.02570097  0.07904955  0.12455686] [0.05555556]
[ 0.5        -0.83333333 -0.16666667] [-0.09323691 -0.02570134  0.07904948  0.12455708] [0.05555556]

We can convert direct lattice quantities, like the atomic position, in crystal coordinates. For instance

atom1 = np.array([0.125,-0.125,-0.125]) # in units of alat
P.ParsersUtils.convert_to_crystal(results.get_lattice(rescale=True),atom1)

array([-0.125, -0.125, -0.125])

The k points can be expressed in crystal coordinates, i.e. in terms of the basis vector of the reciprocal lattice.

For instance we compute the first 10 k points in these coordinates

rec_lat = results.get_reciprocal_lattice(rescale=True)
kpoints_crystal = []
for k in results.kpoints[:10]:
    kpoints_crystal.append(P.ParsersUtils.convert_to_crystal(rec_lat,k))
kpoints_crystal

[array([0., 0., 0.]),
 array([0.        , 0.        , 0.16666667]),
 array([0.        , 0.        , 0.33333333]),
 array([ 0. ,  0. , -0.5]),
 array([0.        , 0.16666667, 0.16666667]),
 array([0.        , 0.16666667, 0.33333333]),
 array([ 0.        ,  0.16666667, -0.5       ]),
 array([ 0.        ,  0.16666667, -0.33333333]),
 array([ 0.        ,  0.16666667, -0.16666667]),
 array([0.        , 0.33333333, 0.33333333])]

We can check in the log of QuantumESPRESSO that these are the k points in crystal coordinates.

We can also perform the inverse procedure, for instance

for k in kpoints_crystal:
    print(P.ParsersUtils.convert_to_cartesian(rec_lat,k))

[0. 0. 0.]
[-0.16666667  0.16666667 -0.16666667]
[-0.33333333  0.33333333 -0.33333333]
[ 0.5 -0.5  0.5]
[0.         0.33333333 0.        ]
[-0.16666667  0.5        -0.16666667]
[ 0.66666667 -0.33333333  0.66666667]
[ 0.5        -0.16666667  0.5       ]
[0.33333333 0.         0.33333333]
[0.         0.66666667 0.        ]

There are also several get methods, for instance

results.get_fermi()

6.117171271888534

The method get_evals return an array with the ks energies for each kpoint. The energies are expressed in eV and the energy of VBM is used as reference. A gap, both direct or indirect, can be set. In this case the energies of the empty ks states is shifted. This procedure does not update the occupation levels of the empty states

results.get_evals(set_gap=1.2)[0]

There are no empty bands. No energy shift has been applied

array([-1.18968067e+01, -1.08435998e-05, -6.32798614e-09,  0.00000000e+00])

results.get_gap()

There are no empty states. Gap cannot be computed.

We test the PwParser using a nscf output file, in this case empty bands are present and the gap can be computed

result_nscf = P.PwParser('QeCalculator_test/bands_8.save/data-file-schema.xml')

Parse file : QeCalculator_test/bands_8.save/data-file-schema.xml

result_nscf.occupations[0]

array([1., 1., 1., 1., 0., 0., 0., 0.])

result_nscf.get_gap()

Indirect gap system
===================
Gap : 0.7360305350920111 eV
Direct gap : 2.5651597869383274 eV

{'gap': 0.7360305350920111,
 'direct_gap': 2.5651597869383274,
 'position_cbm': 12,
 'positon_vbm': 0}

result_nscf.get_evals(set_direct_gap=3.0)[0]

Apply a scissor of 0.43484021306167264 eV

array([-1.18968473e+01, -3.77496830e-06, -2.75945489e-10,  0.00000000e+00,
        3.00000000e+00,  3.00000000e+00,  3.00000201e+00,  3.58149154e+00])

we see that in this case the energy of the 5-th states has been shift to 3 eV to implement the set_direct_gap requirement.

It is also possible to add an explicit scissor to shift the empty bands

result_nscf.get_evals(set_scissor=1.)[0]

Apply a scissor of 1.0 eV

array([-1.18968473e+01, -3.77496830e-06, -2.75945489e-10,  0.00000000e+00,
        3.56515979e+00,  3.56515979e+00,  3.56516180e+00,  4.14665133e+00])

We compute the transition energies from the full to empty bands. Here we show the transition for the first k point

result_nscf.get_transitions(set_direct_gap=3.0)[0]

array([14.89684725, 14.89684725, 14.89684927, 15.47833879,  3.00000377,
        3.00000377,  3.00000579,  3.58149531,  3.        ,  3.        ,
        3.00000201,  3.58149154,  3.        ,  3.        ,  3.00000201,
        3.58149154])

Parse of a Graphene output file

Run the Analysis_BandStructure notebook to produce the xml file used by the PwParser.

We parse a graphene output file to test if systems without a gap are correctly parsed.

results_gra = P.PwParser('Pw_bands/graphene_nscf.save/data-file-schema.xml')

Parse file : Pw_bands/graphene_nscf.save/data-file-schema.xml

results_gra.occupations_kind

'smearing'

results_gra.occupations[0]

array([1.00000000e+00, 1.00000000e+00, 1.00000000e+00, 1.00000000e+00,
       4.22181265e-29, 6.89778699e-37, 1.76721673e-41, 1.28104982e-57])

results_gra.occupations[18] #the position of K

array([1. , 1. , 1. , 0.5, 0.5, 0. , 0. , 0. ])

results_gra.get_gap()

Direct gap system
=================
Gap : 2.715685454290906e-10 eV

{'gap': 2.715685454290906e-10,
 'direct_gap': 2.715685454290906e-10,
 'position_cbm': 18,
 'positon_vbm': 18}