This file is legacy

Use the scatteringCalc.py or debye.py or multiJob_debye.py instead



In [ ]:

    
from __future__ import division
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import sasmol.sasmol as sasmol
from compiledUtils._cLoops import ffi, lib
from lennard_gofr import *
import matplotlib as mpl
from GV import *
import cPickle as pickle

# Style plots
plt.style.use('ggplot')
mpl.rcParams['figure.figsize']=(16,9)
mpl.rcParams['font.size']=20
mpl.rcParams['axes.labelsize']=25
mpl.rcParams['axes.titlesize']=25
mpl.rcParams['figure.titlesize']=25

Calculate g(r)



In [ ]:

    
mol = sasmol.SasMol(0)
mol.read_pdb('/home/data/sascalc_pbc/pdbs/run_0.pdb')
mol.read_dcd('/home/data/sascalc_pbc/pdbs/run_1.dcd')

gc = gofr_calc(mol)


for i in xrange(901,910):
    if i%50==0: print(str(i))
    gc.g_hist(i)
hist = np.copy(gc.g)
r,gr = gc.g_of_r()
plt.xlabel(r'r ($\AA}$)')
plt.ylabel('g(r)')
ax = plt.gca()
# plt.title('g(r) full argon sim')

plt.plot(r,gr)
plt.show()

Calculate Scattering Averaged over frames

Utility Function



In [ ]:

    
def cast_matrix(matrix, ffi):
    ap = ffi.new("double* [%d]" % (matrix.shape[0]))
    ptr = ffi.cast("double *", matrix.ctypes.data)
    for i in range(matrix.shape[0]):
        ap[i] = ptr + i*matrix.shape[1]                                                                
    return ap

Settings



In [ ]:

    
startFrame = 1
endFrame   =  -1 #-1 = use all
NUM_Q      = 10
START_Q    = -1
END_Q      = 1.6
N_GV       = 35
gv = GV(N_GV).get_golden_vectors()
Q_list = np.logspace(START_Q,END_Q,NUM_Q)

import datetime
print(datetime.datetime.now())



In [ ]:

    
coor=mol.coor()[startFrame:endFrame]
num = len(coor[0])

# I = np.zeros((len(coor),len(Q_list)))
for frame in xrange(len(coor)):
    if(frame%50==0):
        print(frame)
    for i,Q in enumerate(Q_list):
        I_tmp = 0
        for g in gv:
            q=g*Q
            cast_coor = cast_matrix(coor[frame],ffi)
            cast_q = ffi.cast('double*',q.ctypes.data)
            I_tmp += lib.sQ(cast_coor,cast_q,num,num)
        I[frame][i] = I_tmp/len(gv)
# pickle.dump( I, open( "multiFrame", "wb" ) )
print(datetime.datetime.now())



In [ ]:

    
plt.loglog(Q_list,I[0])



In [ ]:

    
I



In [ ]:

    
pickle.dump( I, open( "multiFrame", "wb" ) )



In [ ]:

    
a=pickle.load(open('multiFrame10-allFrames','rb'))
a



In [ ]:

    
startFrame = 1
endFrame   =  -1 #-1 = use all
NUM_Q      = 100
START_Q    = -1
END_Q      = 1
N_GV       = 35
gv = GV(N_GV).get_golden_vectors()
Q_list = np.logspace(START_Q,END_Q,NUM_Q)

import datetime
print(datetime.datetime.now())



In [ ]:

    
from multiprocessing import Pool
data_inputs = [0]*4
coor=mol.coor()[startFrame:endFrame]
num = len(coor[0])
N_Q = 100
def process_frame(frame):
    I = np.zeros(len(Q_list))
    for i,Q in enumerate(Q_list):
        I_tmp = 0
        for g in gv:
            q=g*Q
            cast_coor = cast_matrix(coor[frame],ffi)
            cast_q = ffi.cast('double*',q.ctypes.data)
            I_tmp += lib.sQ(cast_coor,cast_q,num,num)
        I[i] = I_tmp/len(gv)
    return I
from functools import partial

partial_frame = partial(process_frame,data_inputs)



In [ ]:

    
# Q_list = np.logspace(START_Q,END_Q,10)
pool = Pool(processes=5)              # process per core
frames = np.arange(600,1000,1)
I_mp=pool.map(process_frame, frames)
pickle.dump( I_mp, open( "multiFrame-"+str(NUM_Q)+"-allFrames", "wb" ) )
print(datetime.datetime.now())



In [ ]:

    
plt.loglog(Q_list,I_mp[0],label='mp')
a=pickle.load(open('/home/data/Outputs/2016-07-16_00-24/multiFrame-2-999_-1-1.6_250','rb'))
Q_list = np.logspace(-1,1.6,250)
plt.loglog(Q_list,a[990],label='Orig')
plt.legend()

Debye no PBC



In [ ]:

    
def pairwise_numpy(X):
    return np.sqrt(((X[:, None, :] - X) ** 2).sum(-1))



In [ ]:

    
def pairwise_numpy(X):
    return np.sqrt(((X[:, None, :] - X) ** 2).sum(-1))
import periodictable.cromermann as ptc
coor = mol.coor()
# pw = pairwise_numpy(coor)
Q_list  = np.logspace(-1,1.6,250)
I  = np.zeros_like(Q_list)
sld  = np.square(ptc.fxrayatq('Ar',Q_list))
for frame in xrange(990,997):
    print(frame)
    pw = pairwise_numpy(coor[frame])
    for i,q in enumerate(Q_list):
        if(i%100==0): print (i)
        I[i] = np.nansum(np.sinc(q*pw))
    plt.loglog(Q_list,I,alpha=.15)



In [ ]:

All Numpy Exp?



In [ ]:

    
coor = mol.coor()[-2]
i, j = np.triu_indices(len(coor), 1)
# i,j=np.indices((len(coor),1))
displacements = coor[i] - coor[j]
print(displacements)



In [ ]:

    
i,j=np.indices((len(coor),1))
len(coor[i]-coor[j])
disp = np.zeros((len(coor),len(coor),3))
for i in xrange(len(coor)):
    for j in xrange(len(coor)):
        disp[i][j] = coor[i]-coor[j]
print()



In [ ]:

    
a=np.array([1,2])
b=np.array([[11,12],[13,14],[13,14]])
# print(np.dot(a,b))
print(' ')
print(np.inner(a,b))
print(1*11+2*12)
print(1*13+2*14)
np.inner([2,4,5],disp).size
print(np.inner([1,1,2],disp)[1])
print(np.inner([1,1,2],disp[1]))
print(np.inner([1,1,2],disp.reshape((2048**2,3))))



In [ ]:

    
tst = disp[0:5,0:5]
# print(tst)
# print(tst.shape)
print(tst.reshape((25,3)))
# print(tst.reshape((25,3)).shape)
q=np.array([2,3,5])
mult = np.inner([[1,1,1],[2,2,2],[3,3,3]],tst.reshape((25,3)))
print(mult.shape)
print(mult)
# print(mult[5])
print(np.sum(np.cos(mult)))
print(np.sum(np.cos(mult[0]))+np.sum(np.cos(mult[1])))



In [ ]:

    
%%time
import numexpr as ne
Q_list = np.logspace(-1,.9,250)
N_GV       = 35
gv = GV(N_GV).get_golden_vectors()
I_exp=np.zeros_like(Q_list)
reshaped = disp.reshape((2048**2,3))*3.4
# print(reshaped)
dotted = np.inner(gv,reshaped)

for i,Q in enumerate(Q_list):
    if(i%10==0): print(i)
    inner = dotted*Q
    I_exp[i] += ne.evaluate("sum(cos(inner))")
#     I_exp[i] += np.sum(np.cos(dotted*Q))
#     for g in gv:
#         q = g*Q
#         inner = np.inner(q,disp.reshape((2048**2,3)))
#         I_exp[i] = ne.evaluate("sum(cos(inner))")
#         inner = np.inner(q,disp)
#         for j in xrange(len(disp)):
#             inner_=inner[j]
            #I_exp[i] = np.sum(ne.evaluate("cos(inner)"))#np.cos(np.inner(q,disp[j])))
#             I_exp[i] += ne.evaluate("sum(cos(inner_))")
I_exp/=len(gv)
print(I_exp)



In [ ]:

    
%%time
coor = mol.coor()[-2]*3.4
num = len(coor)
I = np.zeros(len(Q_list))
for i,Q in enumerate(Q_list):
    if (i%5==0): print(i)
    I_tmp = 0
    for g in gv:
        q=g*Q
        cast_coor = cast_matrix(coor,ffi)
        cast_q = ffi.cast('double*',q.ctypes.data)
        I_tmp += lib.sQ(cast_coor,cast_q,num,num)
    I[i] = I_tmp/len(gv)



In [ ]:

    
plt.loglog(Q_list,I,label='C',lw=5)
plt.loglog(Q_list,I_exp,'k-',label = 'Pure Numpy',lw=2)
plt.legend()
plt.show()



In [ ]:

    
print(disp.shape)



In [ ]:

    
%timeit ne.evaluate("sum(cos(inner))")



In [ ]:

    
%timeit np.sum(np.cos(inner))



In [ ]:

    
%%time
import numexpr as ne
len(disp[0])
Q_list = np.logspace(-1,.5,20)
N_GV       = 35
gv = GV(N_GV).get_golden_vectors()
I_exp_old=np.zeros_like(Q_list)
for i,Q in enumerate(Q_list):
    if(i%5==0): print(i)
    for g in gv:
        q = g*Q
        inner = np.inner(q,disp*3.4)
        for j in xrange(len(disp)):
            inner_=inner[j]
            #I_exp[i] = np.sum(ne.evaluate("cos(inner)"))#np.cos(np.inner(q,disp[j])))
            I_exp[i] += ne.evaluate("sum(cos(inner_))")
I_exp_old/=len(gv)
print(I_exp_old)



In [ ]:

    
import numpy as np
import numexpr as ne
import sasmol.sasmol as sasmol
mol = sasmol.SasMol(0)
mol.read_pdb('Data/run_0.pdb')
mol.read_dcd('Data/run_1.dcd')
from GV import *

coor = mol.coor()[900]
Q_list = np.logspace(-1,1.6,50)
N_GV       = 35
gv = GV(N_GV).get_golden_vectors()

disp = np.zeros((len(coor),len(coor),3))
for i in xrange(len(coor)):
    for j in xrange(len(coor)):
        disp[i][j] = coor[i]-coor[j]



In [ ]:

    
def neTimeLoop():
    Q=1
    I_exp = 0
    for g in gv:
        q = g*Q
        inner = np.inner(q,disp)
        for j in xrange(len(disp)):
            inner_=inner[j]
            #I_exp[i] = np.sum(ne.evaluate("cos(inner)"))#np.cos(np.inner(q,disp[j])))
            I_exp += ne.evaluate("sum(cos(inner_))")
    return I_exp

%timeit neTimeLoop()



In [ ]:

    
def neTimeLoop2():
    Q=1
    I_exp = 0
    for g in gv:
        q = g*Q
        inner = np.inner(q,disp)
        I_exp += ne.evaluate("sum(cos(inner))")
#         inner_=inner[5]
#         for j in xrange(len(disp)):
#             inner_=inner[j]
            #I_exp[i] = np.sum(ne.evaluate("cos(inner)"))#np.cos(np.inner(q,disp[j])))

    return I_exp
%timeit neTimeLoop2()



In [ ]:

    
def npTimeLoop(Q=1):
    I_exp = 0
    for g in gv:
        q = g*Q
        inner = np.inner(q,disp)
        for j in xrange(len(disp)):
            I_exp += np.sum(np.cos(inner[j]))
    return I_exp
%timeit npTimeLoop(1)



In [ ]:

    
def npTimeLoop2(Q=1):
    I_exp = 0
    for g in gv:
        q = g*Q
#         inner = np.inner(q,disp)
        I_exp += np.sum(np.cos(np.inner(q,disp)))
#         for j in xrange(len(disp)):
#             I_exp += np.sum(np.cos(inner[j]))
    return I_exp
%timeit npTimeLoop2(1)



In [ ]:

    
from numba import jit
@jit
def npTimeLoopJit(Q=1):
    I_exp = 0
    for g in gv:
        q = g*Q
#         inner = np.inner(q,disp)
        I_exp += np.sum(np.cos(np.inner(q,disp)))
#         for j in xrange(len(disp)):
#             I_exp += np.sum(np.cos(inner[j]))
    return I_exp
%timeit npTimeLoopJit(1)



In [ ]:

    
from _cLoops import ffi, lib
def cast_matrix(matrix, ffi):
    ap = ffi.new("double* [%d]" % (matrix.shape[0]))
    ptr = ffi.cast("double *", matrix.ctypes.data)
    for i in range(matrix.shape[0]):
        ap[i] = ptr + i*matrix.shape[1]
    return ap
def cffiLoop(Q=1):
    I_tmp = 0
#     frame=900
    num=len(coor)
    for g in gv:
        q=g*Q
#         print(q)
#         print(coor[frame])
        cast_coor = cast_matrix(coor,ffi)
        cast_q = ffi.cast('double*',q.ctypes.data)
        I_tmp += lib.sQ(cast_coor,cast_q,num,num)
    return I_tmp
%timeit cffiLoop()



In [ ]:

    
def npTimeLoop2(Q=1):
    I_exp = 0
    for g in gv:
        q = g*Q
#         inner = np.inner(q,disp)
        I_exp += np.sum(np.cos(np.inner(q,disp)))
#         for j in xrange(len(disp)):
#             I_exp += np.sum(np.cos(inner[j]))
    return I_exp
def npExp(Q):
    I = 0
    for g in gv:
        q = g*Q
        I += np.sum(np.cos(np.inner(q,disp)))
    return I/len(gv)
I_exp = np.zeros_like(Q_list)
for i,Q in enumerate(Q_list):
    if(i%5 ==0): print(i)
    I_exp[i] = npExp(Q)



In [ ]:

    
%matplotlib inline
import matplotlib.pyplot as plt
import matplotlib as mpl
plt.style.use('ggplot')
mpl.rcParams['figure.figsize']=(16,9)
mpl.rcParams['font.size']=20
mpl.rcParams['axes.labelsize']=25
mpl.rcParams['axes.titlesize']=25
mpl.rcParams['figure.titlesize']=25



In [ ]:

    
plt.loglog(Q_list,I_exp,label='numpy Exp')
plt.show()



In [ ]:

    
q=[.3,.1,.1]
inner = np.inner(q,disp)
print(inner)
print('\n\n')
print(np.sum(np.cos(inner)))
print(' ')
print(np.apply_along_axis)



In [ ]:

    
i, j = np.triu_indices(len(coor), 0)
print(i,j)
print(len(i))
print(len(j))



In [ ]:

    
displacements



In [ ]:

    
Q_list = np.logspace(-1,1.6,25)
I_exp =  np.zeros_like(Q_list)
for i,Q in enumerate(Q_list):
    if(i%25==0): print(i)
    for g in gv:
        q = g*Q
        I_exp[i] += np.sum(np.cos(np.inner(q,displacements)))
I_exp /= len(gv)
plt.loglog(Q_list,I_exp)



In [ ]:

    
plt.loglog(Q_list,-I_exp)



In [ ]:

    
%timeit np.sum(np.cos(np.inner(q,displacements)))



In [ ]:

    
np.dot(q,np.transpose(displacements))



In [ ]:

    
Q_list



In [ ]:

    
print(type(np.mat(q)))



In [ ]:

    
q=np.array([1,2,3])
disp  = np.array([[1,1,1],[2,2,2],[3,3,3],[3,3,3]])
for i in displacements:
    print(np.dot(q,i))



In [ ]:

    
np.inner(q,displacements)



In [ ]:

    
disp



In [ ]:

    
displacements



In [ ]:

    
np.dot(q,displacements)



In [ ]:

    
import numexpr as ne



In [ ]: