In [17]:
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 [52]:
class atom(object):
def __init__(self,r,element,q,p):
self.element=element
self.r=r
self.q=q
self.p=p
def add_mapq(self,mapq):
self.mapq=mapq
class molecule(object):
def __init__(self):
self.N=False
self.D=False
self.Q=False
self.Dmap=False
self.atomlist=[]
def readfile(self,filename):
line1=False
line2=False
with open(filename,"r") as f:
lines=f.readlines()
for line in lines:
a=line.split()
if line[0]=="!" and "N=" in line:
self.Q=float(a[2].split("=")[1])
self.D=np.array([float(a[3].split("=")[1]),float(a[4]),float(a[5])])
elif "Rank" in line:
element=a[0]
#print a[1:4]
r=0.1*np.array(a[1:4],dtype=float) #A to nm
line1=True
elif len(a)==1 and line1==True:
q=float(a[0])
line2=True
elif "P" in line and line1==True and line2==True:
p=np.array(a[1:],dtype=float)
self.atomlist.append(atom(r,element,q,p))
def mapcharges(self,molecule):
for true,mapped in zip(self.atomlist,molecule.atomlist):
if true.element==mapped.element:
true.add_mapq(mapped.q)
def calcCofM(self):
cofM=[]
for i in self.atomlist:
cofM.append(i.r)
#print cofM
self.cofM=np.average(cofM,axis=0)
def calcDip(self):
Dip=[]
self.calcCofM()
for i in self.atomlist:
Dip.append(i.q*(i.r-self.cofM))
self.Dip=np.sum(Dip,axis=0)
def calcDipmapped(self,molecule):
self.mapcharges(molecule)
Dip=[]
self.calcCofM()
for i in self.atomlist:
Dip.append(i.mapq*(i.r-self.cofM))
self.Dipmapped=np.sum(Dip,axis=0)
if np.dot(self.Dipmapped,self.Dip)<0:
self.Dipmapped=-self.Dipmapped
In [53]:
#files1=os.listdir("logfiles2")
mapfile="changedorderDCV2T.mps"
mapmol=molecule()
mapmol.readfile(mapfile)
files=[]
for i in range(1,513):
files.append("DCV2T_n2s1_{}.mps".format(i))
#files2=os.listdir("espfits")
molecules=[]
for filename in files:
temp=molecule()
temp.readfile("espfits/"+filename)
temp.calcDip()
temp.calcDipmapped(mapmol)
molecules.append(temp)
trdipcl=[]
trdipmap=[]
for mol in molecules:
trdipcl.append(mol.Dip)
trdipmap.append(mol.Dipmapped)
In [53]:
In [63]:
orarray=np.array(trdipcl)
prarray=np.array(trdipmap)
#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")
for c,q in zip(orarray,prarray):
#if np.absolute(c[5]-q[5])<np.pi and np.absolute(c[4]-q[4])<np.pi/2.0:
x=[c[0],q[0]]
y=[c[1],q[1]]
z=[c[2],q[2]]
#print x,y
ax.plot(x,y,z,c="green")
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 [62]:
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):
#if np.absolute(c[5]-q[5])<np.pi and np.absolute(c[4]-q[4])<np.pi/2.0:
x=[c[5],q[5]]
y=[c[4],q[4]]
#print x,y
plt.plot(x,y,c="green")
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 [74]:
a=plt.hist(cl[:,3], 30, normed=1, facecolor='red', alpha=0.75)
b=plt.hist(qm[:,3], 30, normed=1, facecolor='blue', alpha=0.75)
plt.show()
print np.average(cl[:,3])**4/np.average(qm[:,3])**4
In [ ]: