

In [36]:

    
import numpy as np
from pprint import pprint
from mpi4py import MPI
import lindblad_phasepoints as lb
from scipy.integrate import odeint
import matplotlib.pyplot as plt
from matplotlib import rc

comm = MPI.COMM_WORLD
rank = comm.Get_rank()
size = comm.Get_size()
amp = 40.0
det = 0.0

#############################################################
##Parameters

##Output file path
outfilepath="/home/daneel/workspace/"

##(Un)comment block below for 2 atoms
lattice_size = 2
##Atom positions
X = [1.56886332,-2.19684984]
Y = [2.13749844,-1.48964869]
Z = [0.01768394,1.50295756]
c = np.vstack((X,Y,Z)).T
##Choose the sites for output
i,j = 0,1

##(Un)comment block below for 8 atoms
#lattice_size = 8
##Atom positions
#c = np.array(\
#  [[2.8905099e+00, -6.4307892e-01, -2.2003016e+00], \
#  [-2.7971095e+00, -5.7052033e+00, -1.5733199e+00], \
#  [-1.3179098e+00, -9.9783672e-01, -4.4932801e+00], \
#  [ 2.4362181e+00,  7.2168340e-01, -2.6250514e-01], \
#  [ 1.9754890e+00,  5.7246455e+00, -1.2107655e+00], \
#  [-1.1571209e+00, -3.4153661e+00,  1.2492316e+00], \
#  [-4.8293769e-01, -1.4840459e+00,  1.3405251e-01], \
#  [-3.6379785e-01, -9.0011327e-01,  2.4887775e+00]])
##Choose the sites for output
#i,j = 2,4
#############################################################


l = lattice_size
#Build the atoms
a = np.array([lb.Atom(coords = c[k], index = k) for k in xrange(l)])


#Initiate the parameters in object
p = lb.ParamData(latsize=l, amplitude=amp, detuning=det)
#Initiate the BBGKY system with the parameters 
d = lb.BBGKY_System_Noneqm(p, comm, atoms=a, verbose=False)

#Prepare the times
t0 = 0.0
tfinal = 15.0
nsteps = 3000
times = np.linspace(t0, tfinal, nsteps)

fullsize = 9*l*l + 3*l

#Initial density matrix
rho_init = np.zeros(fullsize)
rho_init[2*l:3*l]=np.ones(l)

#BBGKY time evolution
rho_t = odeint(lb.lindblad_bbgky_pywrap, rho_init, times, args=(d,), Dfun=None)

#Plot gxx_ij
rc('text', usetex=True)
fs=16
gxx_ij_t = rho_t[:,3*l + (j+l*i)]
gyz_ij_t = rho_t[:,3*l + (2 + 3 * 1) * l * l + (j + l * i)]

plt.subplot(2, 1, 1)
plt.plot(times,gxx_ij_t.real,'b-',label="g^{xx}_{ij}")
plt.plot(times,gyz_ij_t.real,'g-',label="g^{yz}_{ij}")
plt.title('(i,j) = '+str((i,j))+', $\Omega$ = '+str(amp)+', $\Delta_0$ = '+str(det)+', N ='+str(l))
plt.legend(loc='best', fancybox=True, framealpha=0.5)
plt.subplot(2, 1, 2)
plt.plot(times, gxx_ij_t.real, 'b')
plt.plot(times,gyz_ij_t.real,'g')
plt.xlim((2.0,2.5))
plt.xlabel('time')

plt.show()

outfilename_xx="gxx_"+str(i)+str(j)+"_time"+"_N_"+str(l)+".txt"
outfilename_yz="gyz_"+str(i)+str(j)+"_time"+"_N_"+str(l)+".txt"
np.savetxt(outfilepath+outfilename_xx,np.vstack((times, gxx_ij_t.real)).T,delimiter=" ")
np.savetxt(outfilepath+outfilename_yz,np.vstack((times, gxx_ij_t.real)).T,delimiter=" ")