In [65]:
import numpy as np
import sys
import matplotlib.pyplot as plt
import os
import scipy.constants as ct
from mpl_toolkits.mplot3d import Axes3D
%matplotlib inline
Ryd2eV=13.605692
bohr2nm=ct.physical_constants["Bohr radius"][0]*1e9
print bohr2nm
plt.rcParams['figure.figsize'] = (14.0, 10.0)
In [51]:
def readmps(filename):
atom1=4
atom2=15
trdip=0
with open(filename,"r") as f:
lines=f.readlines()
if "D[e*nm]" in lines[1]:
a=lines[1].split()
z=float(a[5])
y=float(a[4])
x=float(a[3].split("=")[1])
trdip=np.array([x,y,z])
b=lines[atom1-1].split()
x=float(b[1])
y=float(b[2])
z=float(b[3])
at1=np.array([x,y,z])
c=lines[atom2-1].split()
x=float(c[1])
y=float(c[2])
z=float(c[3])
at2=np.array([x,y,z])
orientation=at1-at2
return trdip,orientation
def readorb(filename):
trdip=0
with open(filename,"r") as f:
lines=f.readlines()
for line in lines:
if "Singlet1 x y z" in line:
#print line
a=line.split()
x=float(a[7])
y=float(a[8])
z=float(a[9])
trdip=-1*np.array([x,y,z])*bohr2nm
break
return trdip
In [52]:
#files1=os.listdir("logfiles2")
files1=[]
files2=[]
for i in range(1,513):
files2.append("DCV2T_n2s1_{}.mps".format(i))
files1.append("molecule_{}.orb.log".format(i))
#files2=os.listdir("espfits")
trdipqm=[]
for filename in files1:
#print filename
trp=readorb("logfiles2/"+filename)
trdipqm.append(trp)
trdipcl=[]
orient=[]
for filename in files2:
trp,orientation=readmps("espfits/"+filename)
orient.append(orientation)
#print orientation
trdipcl.append(trp)
In [53]:
same=0
diff=0
plus=0
minus=0
dotplus=0
dotminus=0
for i,j in zip(trdipqm,trdipcl):
#print i ,j
dot=np.dot(i,j)
if dot>0:
dotplus+=1
else:
dotminus+=1
if (i[1]>0 and j[1]<0):
plus+=1
elif (i[1]<0 and j[1]>0):
minus+=1
if (i[1]>0 and j[1]>0) or (i[1]<0 and j[1]<0):
same+=1
else:
diff+=1
print same, diff
print plus, minus
print dotplus, dotminus
In [66]:
samedir=0
oppodir=0
for i,j in zip(orient,trdipcl):
#print i,j
dot=np.dot(i,j)
if dot>0:
samedir+=1
else:
oppodir+=1
print samedir, oppodir
orarray=np.array(trdipcl)
prarray=np.array(trdipqm)
#print orarray.T[0]
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.scatter(orarray.T[0], orarray.T[1],orarray.T[2], marker="o",c="b")
ax.scatter(prarray.T[0], prarray.T[1],prarray.T[2], marker="x",c="r")
max_range = np.array([orarray.T[0].max()-orarray.T[0].min(), orarray.T[1].max()-orarray.T[1].min(), orarray.T[2].max()-orarray.T[2].min()]).max()
Xb = 0.5*max_range*np.mgrid[-1:2:2,-1:2:2,-1:2:2][0].flatten() + 0.5*(orarray.T[0].max()+orarray.T[0].min())
Yb = 0.5*max_range*np.mgrid[-1:2:2,-1:2:2,-1:2:2][1].flatten() + 0.5*(orarray.T[1].max()+orarray.T[1].min())
Zb = 0.5*max_range*np.mgrid[-1:2:2,-1:2:2,-1:2:2][2].flatten() + 0.5*(orarray.T[2].max()+orarray.T[2].min())
# Comment or uncomment following both lines to test the fake bounding box:
for xb, yb, zb in zip(Xb, Yb, Zb):
ax.plot([xb], [yb], [zb], 'w')
plt.grid()
plt.show()
In [55]:
def appendSpherical_np(xyz):
ptsnew = np.hstack((xyz, np.zeros(xyz.shape)))
xy = xyz[:,0]**2 + xyz[:,1]**2
ptsnew[:,3] = np.sqrt(xy + xyz[:,2]**2)
ptsnew[:,4] = np.arctan2(np.sqrt(xy), xyz[:,2]) # for elevation angle defined from Z-axis down
#ptsnew[:,4] = np.arctan2(xyz[:,2], np.sqrt(xy)) # for elevation angle defined from XY-plane up
ptsnew[:,5] = np.arctan2(xyz[:,1], xyz[:,0])
return ptsnew
In [68]:
cl=appendSpherical_np(orarray)
qm=appendSpherical_np(prarray)
#plt.scatter(cl[:,5],cl[:,4],c="b",marker="o")
#plt.scatter(qm[:,5],qm[:,4],c="r",marker="x")
for c,q in zip(cl,qm):
x=[c[5],q[5]]
y=[c[4],q[4]]
#print x,y
plt.plot(x,y,lw=2)
x=np.linspace(-np.pi,np.pi,9)
labels=x/np.pi*180
plt.xticks(x, labels)
plt.yticks(x[4:], labels[4:])
plt.show()
In [ ]: