I want to plot the excitonic spectrum for 2D black phosphorus, as well as the orbital component for each energy level in spectrum with yambo-py code. I refer to this website
http://www.yambo-code.org/wiki/index.ph ... e_excitons. Firstly, after getting two files named 'o-bse_y.exc_E_sorted' and 'o-bse_y.exc_I_sorted', I run my code based on yambopy as following,
Code: Select all
import sys
from yambopy import *
a = YamboBSEAbsorptionSpectra(job_string='bse_y',path='.')
excitons = a.get_excitons(min_intensity=0.01,max_energy=3.5,Degen_Step=0.01)
print( "nexcitons: %d"%len(excitons) )
print( "excitons:" )
print( " Energy Intensity Index")
for exciton in excitons:
print( "%8.4lf %8.4lf %5d"%tuple(exciton) )
a.get_wavefunctions(Degen_Step=0.01,repx=range(-1,2),repy=range(1),repz=range(-1,2),Cells=[1,1,1],Hole=[0,24.45,7.80], FFTGvecs=10,wf=True, Direction='13')
a.write_json()
Secondly, I use my code to read .json file to plot orbital figures for energy level. But the .xsf file I got shows nothing, only blank. This part I use the code as following,
Code: Select all
import matplotlib
matplotlib.use('Agg') # prevents crashes if no X server present
from yambopy import *
from qepy import *
import json
import matplotlib.pyplot as plt
import numpy as np
from math import sqrt, ceil
import argparse as args
def get_var(dictionary,variables):
for var in variables:
if var in dictionary:
return dictionary[var]
raise ValueError( 'Could not find the variables %s in the output file'%str(variables) )
f = open('absorptionspectra.json')
data = json.load(f)
f.close()
print 'JSON file loaded.'
print 'Absorption spectra ...'
plt.plot(get_var(data,['E/ev','E/ev[1]']), get_var(data,['EPS-Im[2]' ]),label='BSE',lw=3)
nexcitons = len(data['excitons'])
print 'Number of excitons: ', nexcitons
for n,exciton in enumerate(data['excitons']):
plt.axvline(exciton['energy'])
plt.xlabel('$\omega$ (eV)')
plt.ylabel('Intensity (arb. units)')
plt.ylim(0.0,30.0,5.0)
plt.legend()
plt.savefig('excitonE_abs.png', bbox_inches='tight')
print 'excitonE_abs.png'
### Lattice : necessary to determine bounds of exciton weights plot
print 'Reciprocal lattice ...'
lat = np.array(data['lattice'])
rlat = rec_lat(lat)
x,y,z = np.array(rlat )
cut=0.2 ##default value setted by tsli
xmin,ymin,_ = -(x+y)*cut
xmax,ymax,_ = +(x+y)*cut
print 'Excitons weights ...'
nx = int(ceil(sqrt(nexcitons)))
ny = int(ceil(nexcitons*1.0/nx))
print "cols:",nx
print "rows:",ny
cmap = plt.get_cmap("gist_heat_r")
fig = plt.figure(figsize=(nx*4,ny*4))
sorted_excitons = sorted(data['excitons'],key=lambda x: x['energy'])
for n,exciton in enumerate(sorted_excitons):
w = np.array(exciton['weights'])
qpt = np.array(exciton['qpts'])
ax = plt.subplot(ny,nx,n+1)
ax.scatter(qpt[:,0], qpt[:,1], s=20, c=w, marker='H', cmap=cmap, lw=1, label="e: %lf"%exciton['energy'])
plt.xlim([-cut,cut])
plt.ylim([-cut,cut])
plt.gca().yaxis.set_major_locator(plt.NullLocator())
plt.gca().xaxis.set_major_locator(plt.NullLocator())
ax.set_aspect('equal')
plt.savefig('mono_ex.png', bbox_inches='tight', dpi=1500)
print 'mono_ex.png'
print 'Done.'
Best regards,
Tianshu Li,
Jilin University, China