from time import clock
import matplotlib.pyplot as plt
import scipy

from qutip import *
from qutip.piqs import *

from scipy.sparse import load_npz, save_npz

N = 20
ntls = N
nds = num_dicke_states(N)
print("The number of Dicke states is", nds)
[jx, jy, jz] = jspin(N)
jp = jspin(N, "+")
jm = jspin(N, "-")
system = Dicke(N = N)

###The number of Dicke states is 121

a = 1/np.sqrt(2)
b = 1/np.sqrt(2)
css_symmetric = css(N, a, b)
css_antisymmetric = css(N, a,-b)
excited_ = dicke(N, N/2,N/2)
superradiant = dicke(N,N/2,0)
subradiant = dicke(N,j_min(N),-j_min(N))
ground_ = dicke(N,N/2,-N/2)
ghz_ = ghz(N)

# here we set the initial coefficients
gE = 0 # local emission
gD = 0 # local dephasing
gP = 0 # local pumping
gCE = 1 # collective emission
gCD = 0 # collective dephasing
gCP = 0 # collective pumping
w0 = 1  # bare frequency
wi = 0  # coherent drive frequency

# spin hamiltonian
h0 = w0 * jz
hint = wi * jx
h = h0 #+ hint

#set initial conditions for spins by initializing the system and building the Liouvillian matrix
system = Dicke(hamiltonian = h, 
               N = N, 
               emission = gE, 
               pumping = gP, 
               dephasing = gD, 
               collective_emission = gCE, 
               collective_pumping = gCP, 
               collective_dephasing = gCD)
clock_t0 = clock()
lind = system.lindbladian()
liouv = system.liouvillian()
clock_tf = clock()
dt_clock = clock_tf - clock_t0
print("Time (in seconds) to generate the Liouvillian for N = 20 TLSs:", dt_clock)

###Time (in seconds) to generate the Liouvillian for N = 20 TLSs: 0.05499199999999993

## Solution of the dynamics for different initial conditions

# parameters for the time integration of the dynamics
nt = 1001
td0 = np.log(N)/(N*gCE) # delay time is used as a reference
tmax = 10 * td0
t = np.linspace(0, tmax, nt)

# initial states 
rho01 = excited_
rho02 = superradiant
rho03 = css_symmetric
rho04 = css_antisymmetric
rho05 = subradiant
rho06 = ghz_

#Excited
clock_t0 = clock()
result1 = mesolve(liouv, rho01, t, [], e_ops = [jz, jp*jm, jz**2], 
                  options = Options(store_states=True))
rhot1 = result1.states
jz_t1 = result1.expect[0]
jpjm_t1 = result1.expect[1]
jz2_t1 = result1.expect[2]
clock_tf = clock()
dt_clock = clock_tf - clock_t0
print("Elapsed time (in seconds) for this run: ", dt_clock)

#Superradiant
clock_t0 = clock()
result2 = mesolve(liouv, rho02, t, [], e_ops = [jz, jp*jm, jz**2], 
                  options = Options(store_states=True))
rhot2 = result2.states
jz_t2 = result2.expect[0]
jpjm_t2 = result2.expect[1]
jz2_t2 = result2.expect[2]
clock_tf = clock()
dt_clock = clock_tf - clock_t0
print("Elapsed time (in seconds) for this run: ", dt_clock)

#CSS Symmetric
clock_t0 = clock()
result3 = mesolve(liouv, rho03, t, [], e_ops = [jz, jp*jm, jz**2], 
                  options = Options(store_states=True))
rhot3 = result3.states
jz_t3 = result3.expect[0]
jpjm_t3 = result3.expect[1]
jz2_t3 = result3.expect[2]
clock_tf = clock()
dt_clock = clock_tf - clock_t0
print("Elapsed time (in seconds) for this run: ", dt_clock)

#CSS Antisymmetric
clock_t0 = clock()
result4 = mesolve(liouv, rho04, t, [], e_ops = [jz, jp*jm, jz**2], 
                  options = Options(store_states=True))
rhot4 = result4.states
jz_t4 = result4.expect[0]
jpjm_t4 = result4.expect[1]
jz2_t4 = result4.expect[2]
clock_tf = clock()
dt_clock = clock_tf - clock_t0
print("Elapsed time (in seconds) for this run: ", dt_clock)

#Subradiant
clock_t0 = clock()
result5 = mesolve(liouv, rho05, t, [], e_ops = [jz, jp*jm, jz**2], 
                  options = Options(store_states=True))
rhot5 = result5.states
jz_t5 = result5.expect[0]
jpjm_t5 = result5.expect[1]
jz2_t5 = result5.expect[2]
clock_tf = clock()
dt_clock = clock_tf - clock_t0
print("Elapsed time (in seconds) for this run: ", dt_clock)

#GHZ
clock_t0 = clock()
result6 = mesolve(liouv, rho06, t, [], e_ops = [jz, jp*jm, jz**2], 
                  options = Options(store_states=True))
rhot6 = result6.states
jz_t6 = result6.expect[0]
jpjm_t6 = result6.expect[1]
jz2_t6 = result6.expect[2]
clock_tf = clock()
dt_clock = clock_tf - clock_t0
print("Elapsed time (in seconds) for this run: ", dt_clock)

###
##Elapsed time (in seconds) for this run:  1.4128720000000001
##Elapsed time (in seconds) for this run:  1.3064499999999994
##Elapsed time (in seconds) for this run:  1.2567909999999998
##Elapsed time (in seconds) for this run:  1.333939
##Elapsed time (in seconds) for this run:  0.5167439999999992
##Elapsed time (in seconds) for this run:  1.3920459999999997

#####Visualization

jmax = (0.5 * N)
j2max = (0.5 * N + 1) * (0.5 * N)

label_size = 20
label_size2 = 20
label_size3 = 20
label_size4 = 15
plt.rc('text', usetex = True)


fig1 = plt.figure()
plt.plot(t/td0, jz_t1/jmax, '-.', label = r"$|\frac{N}{2},\frac{N}{2}\rangle$")
plt.plot(t/td0, jz_t2/jmax, '-', label = r"$|\frac{N}{2},0\rangle$")
plt.plot(t/td0, jz_t3/jmax, '-', label = r"$|+\rangle_\mathrm{CSS}$")
plt.plot(t/td0, jz_t4/jmax, '-', label = r"$|-\rangle_\mathrm{CSS}$")
plt.plot(t/td0,jz_t5/jmax, '-', label = r"$|0,0\rangle$")
plt.plot(t/td0,jz_t6/jmax, '--', label = r"$|\mathrm{GHZ}\rangle$")
#plt.ylim([-1,1])
plt.xlabel(r'$t/t_\mathrm{D}$', fontsize = label_size3)
plt.ylabel(r'$\langle J_{z} \rangle(t)$', fontsize = label_size3)
plt.xticks([0,(tmax/2)/td0,tmax/td0])
plt.yticks([-1,0,1])
plt.legend(fontsize = label_size4, ncol = 2)

# plot
plt.title(r'Total Inversion', fontsize = label_size3)    
plt.show()
plt.close()

fig2 = plt.figure()
plt.plot(t/td0, jpjm_t1/j2max, '-.', label = r"$|\frac{N}{2},\frac{N}{2}\rangle$")
plt.plot(t/td0, jpjm_t2/j2max, '-', label = r"$|\frac{N}{2},0\rangle$")
plt.plot(t/td0, jpjm_t3/j2max, '-', label = r"$|+\rangle_\mathrm{CSS}$", linewidth = 2)
plt.plot(t/td0, jpjm_t4/j2max, '-', label = r"$|-\rangle_\mathrm{CSS}$")
plt.plot(t/td0,jpjm_t5/j2max, '-', label = r"$|0,0\rangle$")
plt.plot(t/td0,jpjm_t6/j2max, '--', label = r"$|\mathrm{GHZ}\rangle$")
#plt.ylim([0,1])
plt.xticks([0,(tmax/2)/td0,tmax/td0])
plt.yticks([0,0.5,1])
plt.xlabel(r'$t/t_\mathrm{D}$', fontsize = label_size3)
plt.ylabel(r'$\langle J_{+}J_{-}\rangle(t)$', fontsize = label_size3)
plt.legend(fontsize = label_size4, ncol = 2)

# plot
plt.title(r'Light Emission', fontsize = label_size3)    
plt.show()
plt.close()

fig3 = plt.figure()
djz1 = (jz2_t1- jz_t1**2)/jmax**2
djz2 = (jz2_t2- jz_t2**2)/jmax**2
djz3 = (jz2_t3- jz_t3**2)/jmax**2
djz4 = (jz2_t4- jz_t4**2)/jmax**2
djz5 = (jz2_t5- jz_t5**2)/jmax**2
djz6 = (jz2_t6- jz_t6**2)/jmax**2

plt.plot(t/td0, djz1, '-.', label = r"$|\frac{N}{2},\frac{N}{2}\rangle$")
plt.plot(t/td0, djz2, '-', label = r"$|\frac{N}{2},0\rangle$")
plt.plot(t/td0, djz3, '-', label = r"$|+\rangle_\mathrm{CSS}$")
plt.plot(t/td0, djz4, '-', label = r"$|-\rangle_\mathrm{CSS}$")
plt.plot(t/td0, djz5, '-', label = r"$|0,0\rangle$")
plt.plot(t/td0, djz6, '--', label = r"$|\mathrm{GHZ}\rangle$")
#plt.ylim([-1,1])
plt.xticks([0,(tmax/2)/td0,tmax/td0])
plt.yticks([0,0.5,1])
plt.xlabel(r'$t/t_\mathrm{D}$', fontsize = label_size3)
plt.ylabel(r'$\Delta J_{z}^2(t)$', fontsize = label_size3)
plt.rc('xtick', labelsize=label_size) 
plt.rc('ytick', labelsize=label_size)
plt.legend(fontsize = label_size4, ncol = 2)
plt.title(r'Second-moment Collective Correlations', fontsize = label_size3)
plt.show()
plt.close()
 
by

Python Online Compiler

Write, Run & Share Python code online using OneCompiler's Python online compiler for free. It's one of the robust, feature-rich online compilers for python language, supporting both the versions which are Python 3 and Python 2.7. Getting started with the OneCompiler's Python editor is easy and fast. The editor shows sample boilerplate code when you choose language as Python or Python2 and start coding.

Taking inputs (stdin)

OneCompiler's python online editor supports stdin and users can give inputs to programs using the STDIN textbox under the I/O tab. Following is a sample python program which takes name as input and print your name with hello.

import sys
name = sys.stdin.readline()
print("Hello "+ name)

About Python

Python is a very popular general-purpose programming language which was created by Guido van Rossum, and released in 1991. It is very popular for web development and you can build almost anything like mobile apps, web apps, tools, data analytics, machine learning etc. It is designed to be simple and easy like english language. It's is highly productive and efficient making it a very popular language.

Tutorial & Syntax help

Loops

1. If-Else:

When ever you want to perform a set of operations based on a condition IF-ELSE is used.

if conditional-expression
    #code
elif conditional-expression
    #code
else:
    #code

Note:

Indentation is very important in Python, make sure the indentation is followed correctly

2. For:

For loop is used to iterate over arrays(list, tuple, set, dictionary) or strings.

Example:

mylist=("Iphone","Pixel","Samsung")
for i in mylist:
    print(i)

3. While:

While is also used to iterate a set of statements based on a condition. Usually while is preferred when number of iterations are not known in advance.

while condition  
    #code 

Collections

There are four types of collections in Python.

1. List:

List is a collection which is ordered and can be changed. Lists are specified in square brackets.

Example:

mylist=["iPhone","Pixel","Samsung"]
print(mylist)

2. Tuple:

Tuple is a collection which is ordered and can not be changed. Tuples are specified in round brackets.

Example:

myTuple=("iPhone","Pixel","Samsung")
print(myTuple)

Below throws an error if you assign another value to tuple again.

myTuple=("iPhone","Pixel","Samsung")
print(myTuple)
myTuple[1]="onePlus"
print(myTuple)

3. Set:

Set is a collection which is unordered and unindexed. Sets are specified in curly brackets.

Example:

myset = {"iPhone","Pixel","Samsung"}
print(myset)

4. Dictionary:

Dictionary is a collection of key value pairs which is unordered, can be changed, and indexed. They are written in curly brackets with key - value pairs.

Example:

mydict = {
    "brand" :"iPhone",
    "model": "iPhone 11"
}
print(mydict)

Supported Libraries

Following are the libraries supported by OneCompiler's Python compiler

NameDescription
NumPyNumPy python library helps users to work on arrays with ease
SciPySciPy is a scientific computation library which depends on NumPy for convenient and fast N-dimensional array manipulation
SKLearn/Scikit-learnScikit-learn or Scikit-learn is the most useful library for machine learning in Python
PandasPandas is the most efficient Python library for data manipulation and analysis
DOcplexDOcplex is IBM Decision Optimization CPLEX Modeling for Python, is a library composed of Mathematical Programming Modeling and Constraint Programming Modeling