In [97]:

    

import sqlite3
import numpy as np
import scipy.linalg as lg
import scipy.integrate as it
import scipy.optimize as op
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 [155]:

    

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
mini=np.min([np.min(rates12),np.min(rates21)])
print mini

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/mini
    matrix[i,i]-=rate12/mini
    matrix[j,i]=rate21/mini
    matrix[j,j]-=rate21/mini
  
initial=np.zeros(dimension)

initial[600]=1.0

b=matrix[:-1,-1,np.newaxis]
reducedmatrix=matrix[:-1,:-1]-b


    

    




5.76568162428e+14
844.160901478
[[ 10.73884   0.        0.     ]
 [  0.       10.78449   0.     ]
 [  0.        0.       10.39411]]
1568
56151

In [156]:

    

def f(t,y): 
    value=np.dot(reducedmatrix,y)+b[:,0]
    return value
def jac(t,y):
    return reducedmatrix


    

In [157]:

    

r = it.ode(f, jac).set_integrator('Isoda', with_jacobian=True,atol=1e-12,rtol=1e-12)
r.set_initial_value(initial[:-1])

t1 = 10e-12*mini
print t1
dt = 10e-15*mini
solution=[]
time=[]
while r.successful() and r.t < t1:
    r.integrate(r.t+dt)
    solution.append(r.y)
    time.append(r.t)


    

    




8.44160901478e-09

In [169]:

    

timearray=np.array(time)
solutionarray1=np.array(solution).T
print np.shape(timearray),np.shape(solutionarray)
plast=1-np.sum(solutionarray1,axis=0)
print np.sum(solutionarray1,axis=0)


    

    




(1001,) (1568, 1001)
[ 0.78383784  0.83796699  0.86659781 ...,  0.99931502  0.99931518
  0.99931535]

In [170]:

    

print 1-1/float(dimension)
print np.shape(positions)
#0.99931525


    

    




0.999362244898
(1568, 3)

In [171]:

    

solutionarray=np.vstack((solutionarray1,plast))

results=np.vstack((time,solutionarray))


    

In [161]:

    

np.savetxt("poft.txt",results,header="Occupationprobability over time", fmt='%.4e')


    

In [162]:

    

print np.sum(solutionarray,axis=0)
x=np.arange(-3,4)
g=2*np.pi*lg.inv(box)
print g
kx,ky,kz=np.meshgrid(x,x,x)


    

    




[ 1.  1.  1. ...,  1.  1.  1.]
[[ 0.58508976 -0.         -0.        ]
 [ 0.          0.58261311 -0.        ]
 [ 0.          0.          0.60449479]]

In [163]:

    

kindices=np.array([kx.flatten(),ky.flatten(),kz.flatten()])

kvectors=np.dot(kindices.T,g)
print kvectors.shape
print positions.shape
print kindices.T


    

    




(343, 3)
(1568, 3)
[[-3 -3 -3]
 [-3 -3 -2]
 [-3 -3 -1]
 ..., 
 [ 3  3  1]
 [ 3  3  2]
 [ 3  3  3]]

In [164]:

    

V=np.linalg.det(box)
print 1/V
N=float(dimension)
print N


    

    




0.000830721991361
1568.0

In [165]:

    

rm=initial*positions.T
#print np.dot(kvectors,positions.T-rm).T
expkrj=np.exp(-1j*np.dot(kvectors,positions.T-rm))
print solutionarray.shape,expkrj.shape
summand=V/N*np.dot(solutionarray.T,expkrj.T)
summand.shape


    

    




(1568, 1001) (343, 1568)






    
Out[165]:





(1001, 343)

In [166]:

    

# right now only valid for one insertion position
if False:
    rm=initial*positions.T
    rm=rm[:,initial>0]
    M=np.sum(initial)
    expkrm=np.average(np.exp(-1j*np.dot(kvectors,rm)),axis=1)

    pk=summand*expkrm


    

In [167]:

    

print summand[:,16]
print summand[:,13]
print np.sum(summand)


    

    




[-0.12858145+0.08272446j -0.10438655+0.065817j   -0.08665656+0.05426612j
 ..., -0.00038479-0.00017362j -0.00038472-0.00017368j
 -0.00038465-0.00017374j]
[-0.21299312-0.01574008j -0.15539150-0.0121788j  -0.12099083-0.01040613j
 ...,  0.00024801+0.00026689j  0.00024817+0.00026688j
  0.00024832+0.00026686j]
(1334.34823652+6.78471442575e-15j)

In [168]:

    

plt.plot(time,summand[:,16].real)
plt.plot(time,summand[:,10].real)
print time[0]
def f(t,D):
    return np.exp(-D*kvectors[10,0]**2*t)
#plt.xscale("log")

popt, pcov = op.curve_fit(f, timearray[:10], summand[:,10].real[:10])
print popt,pcov
plt.plot(time,f(timearray,popt[0]))
plt.show()


    

    




8.44160901478e-12
[ 1.] inf






    

In [ ]: