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


    

    




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

In [3]:

    




    

In [3]:

    




    

In [4]:

    

def f(t,y):
    return np.dot(matrix,y)
def jac(t,y):
    return matrix
print np.sum(initial)


    

    




1.0

In [173]:

    

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)


    

    




1e-12

In [173]:

    




    

In [174]:

    

timearray=np.array(time)
solutionarray=np.array(solution).T
print np.shape(timearray),np.shape(solutionarray)


    

    




(10001,) (1568, 10001)

In [175]:

    

print np.shape(positions)


    

    




(1568, 3)

In [176]:

    

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


    

In [177]:

    

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


    

In [178]:

    

x=np.arange(-1,2)
g=2*np.pi*lg.inv(box)
print g
kx,ky,kz=np.meshgrid(x,x,x)


    

    




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

In [179]:

    

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

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


    

    




(3, 27)
(1568, 3)
[ 0.  0.  0.]

In [180]:

    

V=np.linalg.det(box)
print 1/V
N=float(dimension)
print N
print np.sum(solutionarray,axis=0)


    

    




0.000830721991361
1568.0
[ 0.96838071  0.93674782  0.90535341 ...,  0.01157011  0.01156992
  0.01156973]

In [181]:

    

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


    

    




[[  0.00000000e+00   0.00000000e+00   0.00000000e+00 ...,   0.00000000e+00
    0.00000000e+00   0.00000000e+00]
 [ -1.27680434e+00  -1.07598560e+00  -8.75166853e-01 ...,   8.75166853e-01
    1.07598560e+00   1.27680434e+00]
 [ -8.63982217e-01  -4.36545646e-01  -9.10907542e-03 ...,   9.10907542e-03
    4.36545646e-01   8.63982217e-01]
 ..., 
 [ -1.58903229e+01  -9.72950624e+00  -3.56868957e+00 ...,   3.56868957e+00
    9.72950624e+00   1.58903229e+01]
 [ -1.67247924e+01  -1.08009685e+01  -4.87714468e+00 ...,   4.87714468e+00
    1.08009685e+01   1.67247924e+01]
 [ -1.74626504e+01  -1.13297451e+01  -5.19683972e+00 ...,   5.19683972e+00
    1.13297451e+01   1.74626504e+01]]






    
Out[181]:





(10001, 27)

In [182]:

    

# 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 [183]:

    

summand[:,13]


    

    
Out[183]:





array([ 0.74343733+0.j,  0.71915238+0.j,  0.69505052+0.j, ...,
        0.00888251+0.j,  0.00888236+0.j,  0.00888222+0.j])

In [168]:

    

plt.plot(time,summand[:,13].real)
plt.xscale("log")


    

In [ ]: