import pylab
pylab.figure(figsize=(7,4))
ax1 = pylab.gca()
ax2 = ax1.twinx()

si_exp_a = 5.4306 # angstrom
si_exp_g = 1.12   # eV
ge_exp_a = 5.6574 # angstrom
ge_exp_g = 0.66   # eV

methods = ['PPS','PBE']
styles  = ['-','--']
labels  = ['Si','SiGe','Ge']

handles = []
for method,style in zip(methods,styles):
    lattice_constants = []
    band_gaps         = []
    for label in labels:
        outfile = "%s_%s.hdf5" % (label,method)
        # Get final lattice constant
        configuration = nlread(outfile, BulkConfiguration)[-1]
        lattice = configuration.bravaisLattice()
        a = lattice.a().inUnitsOf(Angstrom)
        lattice_constants.append(a)
        # Get band gap
        bandstructure = nlread(outfile, Bandstructure)[-1]
        gap = bandstructure.indirectBandGap().inUnitsOf(eV)
        band_gaps.append(gap)
        # Plot results
    x = range(len(band_gaps))
    h1, = ax1.plot(x, band_gaps, 'ro%s' % style, ms=8)
    h2, = ax2.plot(x, lattice_constants, 'bs%s' % style, ms=8)
    handles.append(h1)
    handles.append(h2)

# Plot experimental data
ax1.plot(0, si_exp_g, 'o', color='k')
ax1.plot(2, ge_exp_g, 'o', color='k')
ax2.plot(0, si_exp_a, 's', color='0.6')
ax2.plot(2, ge_exp_a, 's', color='0.6')
ax1.annotate('black dots    : Exp. gaps', xy=(0.3,0.95), xycoords='axes fraction', fontsize='small')
ax1.annotate('grey squares: Exp. lattice constants', xy=(0.3,0.9), xycoords='axes fraction', fontsize='small')

# Legend
txt = ['PPS gaps', 'PPS lattice', 'PBE gaps', 'PBE lattice']
#       'Si exp. gap', 'Ge exp. gap', 'Si exp. l.c.', 'Ge exp. l.c.']
leg = ax1.legend(handles,txt,ncol=2,loc='lower center',fontsize='small')

# Plot details
ax1.set_xticks(x)
ax1.set_xticklabels(labels)
ax1.set_ylabel('Indirect band gap (eV)')
ax2.set_ylabel('Lattice constant (Angstrom)')
pylab.xlim((-0.1, 2.1))
pylab.tight_layout()
pylab.savefig('SiGe.png')
pylab.show()