In [2]:
import sqlite3
import numpy as np
import scipy.linalg as lg
import scipy.integrate as it
import matplotlib.pyplot as plt
import scipy.constants as ct
import copy
from matplotlib.mlab import griddata
from matplotlib.colors import LogNorm
from matplotlib import ticker
%matplotlib inline
plt.rcParams['figure.figsize'] = (10.0, 8.0)
hbar=ct.physical_constants["Planck constant over 2 pi in eV s"][0]
T=300
kbT=T*ct.physical_constants["Boltzmann constant in eV/K"][0]
In [3]:
sqlname="system_individualmps.sql"
sqlstatement="SELECT pairs.seg1, pairs.seg2, pairs.Jeff2s, seg1.eSinglet, seg1.UnXnNs, seg1.UxNxXs, seg2.eSinglet, seg2.UnXnNs, seg2.UxNxXs FROM pairs JOIN segments seg1 ON seg1._id =pairs.seg1 JOIN segments seg2 ON seg2._id =pairs.seg2"
con = sqlite3.connect(sqlname)
with con:
cur = con.cursor()
cur.execute(sqlstatement)
rows = cur.fetchall()
sql=np.array(rows)
lowerlimit=0
reorg12=sql[:,4]+sql[:,8]
dG12=-sql[:,3]+sql[:,6]
reorg21=sql[:,5]+sql[:,7]
dG21=-dG12
rates12=2*np.pi/hbar*sql[:,2]/np.sqrt(4*np.pi*kbT*reorg12)*np.exp(-(dG12+reorg12)**2/(4*reorg12*kbT))
rates21=2*np.pi/hbar*sql[:,2]/np.sqrt(4*np.pi*kbT*reorg21)*np.exp(-(dG21+reorg21)**2/(4*reorg21*kbT))
maxi=np.max([np.max(rates12),np.max(rates21)])
print maxi
sqlstatement="SELECT box11,box12,box13,box21,box22,box23,box31,box32,box33 from frames"
con = sqlite3.connect(sqlname)
with con:
cur = con.cursor()
cur.execute(sqlstatement)
vecs = cur.fetchall()
box=np.array(vecs).reshape((3,3))
print box
sqlstatement="SELECT posX,posY,posZ from segments"
con = sqlite3.connect(sqlname)
with con:
cur = con.cursor()
cur.execute(sqlstatement)
rows2 = cur.fetchall()
positions=np.array(rows2)
#dimension=int(max(np.max(sql[:,0]),np.max(sql[:,1])))
dimension=1568
print dimension
matrix=np.zeros((dimension,dimension))
print len(rows)
for k in range(len(rows)):
rate12=rates12[k]
rate21=rates21[k]
row=rows[k]
i=row[0]-1
j=row[1]-1
matrix[i,j]=rate12
matrix[i,i]-=rate12
matrix[j,i]=rate21
matrix[j,j]-=rate21
initial=np.zeros(dimension)
initial[0]=1.0
In [4]:
if False:
with open("matrix.txt",'w') as f:
for i in range(1,dimension+1):
for j in range(1,dimension+1):
if matrix[i-1,j-1]!=0.0:
f.write("{:6d} {:6d} {:+1.17e}\n".format(i,j,matrix[i-1,j-1]))
In [5]:
if False:
sqlname="system_individualmps.sql"
sqlstatement = "SELECT pairs.id,pairs.seg1, pairs.seg2, pairs.rate12s, pairs.rate21s FROM pairs"
con = sqlite3.connect(sqlname)
with con:
cur = con.cursor()
cur.execute(sqlstatement)
rows = cur.fetchall()
sql=np.array(rows)
#print rows[0]
print max(np.max(sql[:,3]),np.max(sql[:,4]))
minimum=min(np.min(sql[:,3]),np.min(sql[:,4]))
print minimum
maximum=max(np.max(sql[:,3]),np.max(sql[:,4]))
dimension=int(max(np.max(sql[:,1]),np.max(sql[:,2])))
print dimension
lowerlimit=0
#print rows[0]
matrix=np.zeros((dimension,dimension))
for row in rows:
i=row[1]-1
j=row[2]-1
if row[3]>lowerlimit:
matrix[i,j]=row[3]
matrix[i,i]-=(row[3])
if row[4]>lowerlimit:
matrix[j,i]=row[4]
matrix[j,j]-=(row[4])
print np.sum(matrix==0)
print (np.sum(matrix!=0)-dimension)/2
In [6]:
def f(t,y):
return np.dot(matrix,y)
def jac(t,y):
return matrix
print np.sum(initial)
In [26]:
r = it.ode(f, jac).set_integrator('Isoda', with_jacobian=True)
r.set_initial_value(initial)
t1 = 10e-13
print t1
dt = 10e-17
solution=[]
time=[]
while r.successful() and r.t < t1:
r.integrate(r.t+dt)
solution.append(r.y)
time.append(r.t)
In [8]:
seg1=[]
seg2=[]
for sol in solution:
seg1.append(sol[0])
seg2.append(sol[700])
#print seg1
plt.plot(time,seg1,c='r')
plt.plot(time,seg2,c='b')
#plt.xscale('log')
plt.show()
In [9]:
timearray=np.array(time)
solutionarray=np.array(solution).T
print np.shape(timearray),np.shape(solutionarray)
In [10]:
print np.shape(positions)
In [11]:
results=np.vstack((time,solutionarray))
In [12]:
np.savetxt("poft.txt",results,header="Occupationprobability over time", fmt='%.4e')
In [13]:
class segment(object):
def __init__(self,id,pos):
self.id=id
self.pos=pos
self.ghost=False
def addpoft(self,poft,time):
self.poft=poft
self.time=time
def shift(self,pos,Ghost=True):
self.pos+=pos
self.ghost=Ghost
In [14]:
def plotcontour(x,y,z,box,seg,time):
x=np.array(x)
y=np.array(y)
z=np.array(z)
numcols, numrows = 100, 100
xi = np.linspace(x.min(), x.max(), numcols)
yi = np.linspace(y.min(), y.max(), numrows)
xi, yi = np.meshgrid(xi, yi)
zi = griddata(x, y, z, xi, yi)
fig, ax = plt.subplots()
im = ax.contourf(xi, yi, zi,locator=ticker.LogLocator(),vmin=1e-8, vmax=1e-2)
ax.scatter(x, y, c=z, s=10,vmin=zi.min(), vmax=zi.max())
fig.colorbar(im)
plt.title("{:1.3e} s".format(time))
# drawing simbox
#print seg.pos
x1=seg.pos[0]-0.5*box[0,0]
y1=seg.pos[1]-0.5*box[1,1]
x2=x1+box[0,0]
y2=y1+box[1,1]
ax.plot([x1, x2], [y1, y1], color='k', linestyle='-', linewidth=2)
ax.plot([x1, x1], [y1, y2], color='k', linestyle='-', linewidth=2)
ax.plot([x2, x2], [y1, y2], color='k', linestyle='-', linewidth=2)
ax.plot([x1, x2], [y2, y2], color='k', linestyle='-', linewidth=2)
return fig
In [15]:
segments =[]
i=0
for row in rows2:
segments.append(segment(i,np.array(row)))
i+=1
for j,i in zip(segments,range(np.shape(solutionarray)[0])):
j.addpoft(solutionarray[i],timearray)
In [16]:
for segment in segments[1:10]:
plt.plot(segment.time,segment.poft)
In [17]:
print box
seg_pbc=[]
a=range(-1,2)
for segment in segments:
for i in a:
for j in a:
for k in a:
if i==0 and j==0 and k==0:
seg_pbc.append(segment)
else:
temp=copy.deepcopy(segment)
vector=i*box[:,0]+j*box[:,1]+k*box[:,2]
temp.shift(vector)
seg_pbc.append(temp)
In [18]:
for i in range(len(seg_pbc)):
if seg_pbc[i].poft[0]>0.5 and seg_pbc[i].ghost==False:
print i,seg_pbc[i].pos
In [19]:
xyplane=[]
for segment in seg_pbc:
if np.abs(segment.pos[2]-seg_pbc[300].pos[2])<0.3:
xyplane.append(segment)
In [25]:
print len(xyplane)
x=[]
y=[]
for seg in xyplane:
x.append(seg.pos[0])
y.append(seg.pos[1])
times=seg_pbc[13].time
j=0
for i in range(0,len(times)/15,1):
#print times[i]
p0=[]
for seg in xyplane:
p0.append(seg.poft[i])
fig=plotcontour(x,y,p0,box,seg_pbc[13],times[i])
plt.savefig("contour_{:04d}_s.png".format(j))
plt.close()
j+=1
In [ ]: