In [107]:
import numpy as np
import lxml.etree as lxml
import matplotlib.pyplot as plt
%matplotlib inline
ryd2ev=13.60569253
plt.rcParams['figure.figsize'] = (10.0, 8.0)
In [108]:
def readexcitonlogfile(filename,dft=False,qp=False,singlets=False,triplets=False):
dftlist=[]
gwa=[]
qp=[]
s=[]
t=[]
homo=None
lumo=None
tbool=False
sbool=False
qpbool=False
dftenergy=None
with open(filename,"r") as f:
for line in f.readlines():
if "QM energy[eV]:" in line and dftenergy==None:
dftenergy=float(line.split()[-1])
elif dft and "S-C" in line and "S-X" in line:
entries=line.split()
dftlist.append(ryd2ev*float(entries[7]))
gwa.append(ryd2ev*float(entries[-1]))
if "HOMO"==entries[2] and homo==None:
homo=int(entries[4])-1
if "LUMO"==entries[2] and lumo==None:
lumo=int(entries[4])-1
elif "====== Diagonalized quasiparticle energies (Rydberg) ======" in line:
qpbool=True
elif qpbool and qp and "PQP" in line and "DQP" in line:
qp.append(ryd2ev*float(line.split()[-1]))
elif "====== triplet energies (eV) ======" in line:
tbool=True
elif tbool and triplets and "T =" in line:
t.append(float(line.split()[7]))
elif "====== singlet energies (eV) ======" in line:
sbool=True
elif sbool and singlets and "S =" in line:
s.append(float(line.split()[7]))
results=[]
if dft:
results.append(np.array([homo,lumo,dftenergy]))
results.append(np.array(dftlist))
results.append(np.array(gwa))
else:
results.append(None)
results.append(None)
results.append(None)
if qp:
results.append(np.array(qp))
else:
results.append(None)
if singlets:
results.append(np.array(s))
else:
results.append(None)
if triplets:
results.append(np.array(t))
else:
results.append(None)
#print results
return results
In [109]:
def getcouplingfromsplit(filename,states):
dft=False
singlets=False
triplets=False
coupling=[]
for state in states:
if "e" in state or "h" in state:
dft=True
if "s" in state:
singlets=True
if "t" in state:
triplets=True
results=readexcitonlogfile(filename,dft=dft,qp=False,singlets=singlets,triplets=triplets)
#print results
for state in states:
stateindex=None
resultsindex=None
if state=="e_dft" or state=="e":
stateindex=results[0][1]
resultsindex=1
elif state=="h_dft" or state=="h":
stateindex=results[0][0]-1
resultsindex=1
elif state=="e_gw":
stateindex=results[0][1]
resultsindex=2
elif state=="h_gw":
stateindex=results[0][0]-1
resultsindex=2
elif "s" in state:
stateindex=2*int(state[1:])-2
resultsindex=4
elif "t" in state:
stateindex=2*int(state[1:])-2
resultsindex=5
else:
print "state {} not known".format(state)
splitt=results[resultsindex][stateindex+1]-results[resultsindex][stateindex]
#if state=="e":
#print results[resultsindex][stateindex+1]/ryd2ev
#print results[resultsindex][stateindex]/ryd2ev
#print 0.5*splitt
coupling.append(0.5*splitt)
return coupling
In [110]:
def readcouplingxml(filename):
parser=lxml.XMLParser(remove_comments=True)
tree=lxml.parse(filename,parser)
root=tree.getroot()
Je=[]
Jh=[]
for pair in root:
homoA=int(pair.get("homoA"))
homoB=int(pair.get("homoB"))
for overlap in pair:
orbA=int(overlap.get("orbA"))
orbB=int(overlap.get("orbB"))
if orbA==homoA and orbB==homoB:
Je.append(np.absolute(float(overlap.text)))
elif orbA==homoA+1 and orbB==homoB+1:
Jh.append(np.absolute(float(overlap.text)))
return [Je,Jh]
In [111]:
def readexcitoncouplingxml(filename,states):
parser=lxml.XMLParser(remove_comments=True)
tree=lxml.parse(filename,parser)
root=tree.getroot()
resultlist=[]
for pair in root:
types=pair[0]
couplings=[]
for state in states:
results=None
if state[0]=="s":
results=types.find("singlets")
elif state[0]=="t":
results=types.find("triplets")
else:
print "state not known"
number=int(state[1:])
#print number
for coupling in results:
noA=int(coupling.get("excitonA"))
noB=int(coupling.get("excitonB"))
if noA+1==number and noB+1==number:
couplings.append(np.absolute(float(coupling.text)))
resultlist.append(couplings)
return resultlist
In [112]:
def readexcitoncoulingclassical(filename):
parser=lxml.XMLParser(remove_comments=True)
tree=lxml.parse(filename,parser)
root=tree.getroot()
results=[]
for pair in root:
Coupling=pair[0]
results.append(float(Coupling.get("jABstatic")))
return results
In [113]:
distances=[2.5,2.6,2.7,2.8,2.9,3,3.1,3.2,3.3,3.4,3.5,3.6,3.7,3.8,3.9,4,4.1,4.2,4.3,4.4,4.5,5,6,7]
In [114]:
Jh_splitt=[]
Je_splitt=[]
Js_splitt=[]
Jt_splitt=[]
Jh_frag=[]
Je_frag=[]
Js_frag=[]
Jt_frag=[]
Js2_classical=[]
for distance in distances:
foldername="{0}_0.00_0.00_{0:1.2f}".format(distance)
print foldername
temp=getcouplingfromsplit("PBE/"+foldername+"/exciton.log",["e","h","s1","t1"])
Jh_splitt.append(temp[1])
Je_splitt.append(temp[0])
Js_splitt.append(temp[2])
Jt_splitt.append(temp[3])
temp=readcouplingxml("PBE/"+foldername+"/coupling.out.xml")
Je_frag.append(temp[0])
Jh_frag.append(temp[1])
temp=readexcitoncouplingxml("PBE/"+foldername+"/excitoncoupling.out.xml",["s1","t1"])
Js_frag.append(temp[0][0])
Jt_frag.append(temp[0][1])
Js2_classical.append(readexcitoncoulingclassical("PBE/"+foldername+"/excitoncoupling_classical.out.xml")[0])
In [115]:
plt.plot(distances,Jh_splitt,marker="o",label="Jh_split")
plt.plot(distances,Jh_frag,marker="o",label="Jh_frag")
plt.plot(distances,Je_splitt,marker="o",label="Je_split")
plt.plot(distances,Je_frag,marker="o",label="Je_frag")
plt.yscale('log')
plt.xlabel("dimer separation [A]")
plt.ylabel("|J| in eV")
plt.legend()
plt.savefig("charge_coupling_vs_distance.png")
plt.show()
In [116]:
plt.plot(distances,Js_splitt,marker="o",label="Js_split")
plt.plot(distances,Js_frag,marker="o",label="Js_frag")
plt.plot(distances,Jt_splitt,marker="o",label="Jt_split")
plt.plot(distances,Jt_frag,marker="o",label="Jt_frag")
#plt.plot(distances,Js2_classical,marker="o",label="Js_classical")
plt.xlabel("dimer separation [A]")
plt.ylabel("|J| in eV")
plt.yscale('log')
plt.legend()
plt.savefig("exciton_coupling_vs_distance.png")
plt.show()
print "Hallo", Jt_splitt[0],len(Jt_frag)
In [117]:
Es=[]
Et=[]
Edft=[]
for distance in distances:
foldername="{0}_0.00_0.00_{0:1.2f}".format(distance)
results=readexcitonlogfile("PBE/"+foldername+"/exciton.log",dft=True,qp=False,singlets=True,triplets=True)
Es.append(results[4]+results[0][2])
Et.append(results[5])#+results[0][2])
Edft.append(results[0][2])
Et=np.array(Et)
Es=np.array(Es)
print Et.shape
print Es.shape
In [125]:
+1
for i,line in enumerate(Et.T[:8]):
plt.plot(distances,line,marker="o",label="T {}".format(i+1))
plt.xlabel("dimer separation [A]")
plt.ylabel("Omega T in eV")
plt.legend()
plt.show()
In [124]:
for i,line in enumerate(Es.T[:10]):
plt.plot(distances,line,marker="o",label="S {}".format(i+1))
plt.xlabel("dimer separation [A]")
plt.ylabel("Omega S in eV")
plt.legend()
plt.show()
In [120]:
plt.plot(distances,Edft)
Out[120]:
In [120]:
In [120]:
In [120]: