PYTHON 5
###### Computation.py Guest on 20th November 2020 05:59:18 PM
1. # -*- coding: UTF-8 -*-
2.
3. import numpy as np
4. np.set_printoptions(threshold=np.nan)
5. import math, sys, time, cmath, os
6. import multiprocessing as mp
7. from spin_states import *
8. #from superoperator import *
9. try:
10.     import qutip as qt
11.
12. except:
13.     print " Some modules are missing, try installing them using \"installModules.py\""
14.     sys.exit()
15. THRUSTY = 0
16. NPROCS = 1
17. UPARROW = u"\u2191"
18. DOWNARROW = u"\u2193"
19.
20.
21. def HST(x):
22.     if x==0:
23.         return 0.5
24.     else:
25.         return int(x>0)
26.
27.
28. def phi(N,I,M,b,gamma, alfa):
29.     return cmath.phase( 2*M + 1 + b - N*gamma + 2* alfa * 1j * ( I*(I+1) - M*(M+1) )**0.5 )
30.
31. def x(M,b):
32.     return -1-(2* M+1)*b
33.
34. def y(N,I,M,b,gamma, alfa):
35.     if alfa == 0:
36.         return 2*M + 1 + b - N*gamma
37.
38.     return ( (2*M + 1 + b - N*gamma)**2 + 4 * alfa**2 * (I-M) * (I+1+M) )**0.5
40.     x = x.real
41.     threshold = 14
42.     if x == 0:
44.
46.         if x > 0:
47.             return 0
48.         elif x < 0:
49.             return -x
50.
52.         return -x
54.         return 0
56.
57. class Cluster(object):
58.     def __init__(self, p):
59.         self.parameters = {"N" : p[0], "gamma" : p[1], "b" : p[2], "alfa" : p[3], "beta" : p[4],"eta" : p[5], "tau_dot" : p[6], "tau_lead" : p[7]}
60.         for p, v in self.parameters.iteritems():
61.             setattr(self, p, v)
62.         self.define_operators()
63.         self.local_hamiltonian = self.get_local_hamiltonian()
64.
65.     def get_ancilla_states(self):
66.         N = self.N
67.         i = int(N)
68.         while i >= 0:
69.             state = {"I" : N/2.,"M" : - i + N/2.,"p" : 1,"state" : qt.basis(N + 1, i)}
70.
71.             yield state
72.             i -= 1
73.     def define_operators(self):
74.         N = self.N
75.         self.identity_central = qt.qeye(2)
76.         self.Sm = qt.sigmam()
77.         self.Sp = qt.sigmap()
78.         self.Sx = 0.5 * qt.sigmax()
79.         self.Sy = 0.5 * qt.sigmay()
80.         self.Sz = 0.5 * qt.sigmaz()
81.
82.         self.identity_ancilla = qt.qeye(N + 1)
83.         self.Ix = qt.jmat((N/2.), "x")
84.         self.Iy = qt.jmat((N/2.), "y")
85.         self.Iz = qt.jmat((N/2.), "z")
86.         self.Isq = (self.Ix**2 + self.Iy**2 + self.Iz**2).tidyup()
87.         self.Ip = qt.jmat((N/2.), "+")
88.         self.Im = qt.jmat((N/2.), "-")
89.
90.         self.S2 = qt.tensor(self.Sx,  self.identity_ancilla) **2 + qt.tensor(self.Sy,  self.identity_ancilla) **2 + qt.tensor(self.Sz,  self.identity_ancilla) **2
91.         self.I2 = qt.tensor(self.identity_central, self.Ix) **2 + qt.tensor(self.identity_central, self.Iy) **2 + qt.tensor(self.identity_central, self.Iz) **2
92.         self.m_dm = lambda m : qt.ket2dm(qt.basis(N + 1, m))
93.         self.s_dm = lambda s : qt.ket2dm(qt.basis(2, s))
94.         self.identity = qt.tensor(self.identity_central, self.identity_ancilla)
95.         self.M_stg = qt.tensor(self.Sz, self.identity_ancilla)- 1/(N) * qt.tensor(self.identity_central, self.Iz)
96.
97.
98.     def get_local_hamiltonian(self):
99.         N, b, gamma, alfa, beta = self.N, self.b, self.gamma, self.alfa, self.beta
100.
101.         H_elec =  - N * gamma * self.Sz
102.         H_nuclei = -b * self.Iz
103.
104.         H_hf = beta * (2 * qt.tensor(self.Sz, self.Iz) + alfa * qt.tensor(self.Sp,  self.Im) + alfa * qt.tensor(self.Sm, self.Ip ) )
105.
106.         H = qt.Qobj((qt.tensor(H_elec, self.identity_ancilla)  + qt.tensor(self.identity_central, H_nuclei) + H_hf ).full())
107.         return H
108.
109.     def get_transverse_hamiltonian(self, Delta):
110.         H_elec =  Delta * self.Sx
111.
112.         H = qt.Qobj((qt.tensor(H_elec, self.identity_ancilla)).full())
113.         return H
114.
115.
116.     def calculate_eigenstates(self, N, b, gamma, alfa, ham):
117.         def get_eigenstate(ps,I,M,N, b,gamma, alfa,p, ham,ns1,ns2 = False):
118.             if ps == "u":
119.                 state = qt.tensor(qt.basis(2, 0), ns1["state"])
120.                 te = -0.5*N*gamma-I*(b-1)
121.                 rep = "|%s >|%g %g %d>"%(UPARROW, ns1["M"], ns1["I"], ns1["p"])
122.             elif ps == "d":
123.                 state = qt.tensor(qt.basis(2, 1), ns1["state"])
124.                 te = 0.5*N*gamma+I*(b+1)
125.                 rep = "|%s >|%g %g %d>"%(DOWNARROW, ns1["M"], ns1["I"], ns1["p"])
126.             elif ps == "+":
127.                 f = phi(N,I,M,b,gamma, alfa)
128.                 state = math.cos(f/2.) * qt.tensor(qt.basis(2, 0),ns1["state"]) + math.sin(f/2.) * qt.tensor(qt.basis(2, 1),ns2["state"])
129.                 te = 0.5* x(M,b)  + 0.5*  y(N,I,M,b,gamma, alfa)
130.                 rep = "%+.3f|%s >|%g %g %d> %+.3f|%s >|%g %g %d>"%(math.cos(f/2.), UPARROW, ns1["M"], ns1["I"], ns1["p"], math.sin(f/2.), DOWNARROW, ns2["M"], ns2["I"], ns2["p"])
131.             elif  ps == "-":
132.                 f = phi(N,I,M,b,gamma, alfa)
133.                 state = -math.sin(f/2.) * qt.tensor(qt.basis(2, 0),ns1["state"]) + math.cos(f/2.) * qt.tensor(qt.basis(2, 1),ns2["state"])
134.                 te = 0.5* x(M,b)  - 0.5* y(N,I,M,b,gamma, alfa)
135.                 rep = "%+.3f|%s >|%g %g %d> %+.3f|%s >|%g %g %d>"%(-math.sin(f/2.), UPARROW, ns1["M"], ns1["I"], ns1["p"], math.cos(f/2.), DOWNARROW, ns2["M"], ns2["I"], ns2["p"])
136.             return {"ps" : ps,"I" : I,"M" : M,"p" : p,"state" : qt.Qobj(state.full()).tidyup().unit(), "theorEnergy": te,"energy" : ham.matrix_element(state.dag(),state).real, "representation" : rep}
137.
138.         eigenstates = []
139.         for ns in self.get_ancilla_states(N):
140.             I, M, p = ns["I"], ns["M"], ns["p"]
141.
142.             if round(M,2) == round(I,2):
143.                 eigenstates.append(get_eigenstate("u",I,M,N,b, gamma, alfa,p, ham,ns))
144.
145.             else:
146.                 for ps in ("-", "+"):
147.                     states2 = [elem for elem in self.get_ancilla_states(N) if (elem["I"] == I and elem["p"] == p and elem["M"] == M + 1)]
148.                     if states2:
149.                         ns2 = states2[0]
150.                     else:
151.                         print "{} {}".format(ns["I"], ns["M"])
152.                     eigenstates.append(get_eigenstate(ps,I,M,N,b, gamma, alfa,p, ham,ns,ns2))
153.             if M == -I:
154.                 eigenstates.append(get_eigenstate("d",I,M,N,b, gamma, alfa,p, ham,ns))
155.
156.         eigenstates = sorted([s for s in eigenstates], key = lambda s : s["energy"])
157.         return     eigenstates
158.
159.
160.
161.
162.
163. class Computation(object):
164.     def __init__(self, parameters):
165.         self.parameters = parameters
166.
167.         self.folder = "try4"
168.         self.cluster_parameters = []
169.         self.clusters = [Cluster(self.parameters.cluster_parameters[c]) for c in xrange(self.parameters.clusters)]
170.         self.identity = qt.tensor([c.identity for c in self.clusters])
171.         self.local_H = self.get_local_hamiltonian()
172.         self.transverse_H = self.get_transverse_hamiltonian()
173.         self.ising_H = self.get_ising_hamiltonian()
174.
175.     def get_total_hamiltonian(self, t):
176.         v = self.parameters.annealing_v
177.         H =  max(0, 1- v*t) * self.transverse_H + min(1, v*t) * self.local_H + min(1, v*t) * self.ising_H
178.         return H
179.
180.     def get_local_hamiltonian(self):
181.         for i, c in enumerate(self.clusters):
182.             vector = [g.local_hamiltonian if i==j else g.identity for j, g in enumerate(self.clusters) ]
183.             if i==0:
184.                 H = qt.Qobj(qt.tensor(vector).full())
185.             else:
186.                 H += qt.Qobj(qt.tensor(vector).full())
187.
188.         return H
189.
190.     def get_operators_compute(self):
191.         ops = []
192.         for i, c in enumerate(self.clusters):
193.             vector = [qt.tensor(g.Sz, g.identity_ancilla) if i==j else g.identity for j, g in enumerate(self.clusters) ]
194.             ops.append(qt.Qobj(qt.tensor(vector).full()) )
195.
196.         return ops
197.
198.     def get_transverse_hamiltonian(self):
199.         for i, c in enumerate(self.clusters):
200.             vector = [g.get_transverse_hamiltonian(self.parameters.Delta)  if i==j else g.identity for j, g in enumerate(self.clusters)]
201.             if i==0:
202.                 H = qt.Qobj(qt.tensor(vector).full())
203.             else:
204.                 H += qt.Qobj(qt.tensor(vector).full())
205.
206.         return H
207.
208.     def get_ising_hamiltonian(self):
209.
210.         for i, c in enumerate(self.clusters):
211.
212.             for j, g in enumerate(self.clusters):
213.                 if i<j:
214.                     vector = [qt.tensor(d.Sz, d.identity_ancilla) if k in (i, j) else d.identity for k, d in enumerate(self.clusters) ]
215.                     if i==0 and j==1:
216.                         H = self.parameters.J[i,j] * qt.Qobj(qt.tensor(vector).full())
217.                     else:
218.                         H += self.parameters.J[i,j] * qt.Qobj(qt.tensor(vector).full())
219.
220.         return H
221.
222.
223.     def get_dm(self, hamiltonian):
224.         states = sorted(zip(hamiltonian.eigenstates()[1], hamiltonian.eigenstates()[0]), key= lambda x: x[1])
225.         gs_energy = states[0][1]
226.         tau_dot = 0
227.         if tau_dot > 0:
228.             dm = (-h_init/tau_dot).expm().matrix_element(eigenstates[0]["state"].dag(), eigenstates[0]["state"]) * qt.ket2dm(eigenstates[0]["state"])
229.             for es in eigenstates[1:]:
230.                 dm += (-h_init/tau_dot).expm().matrix_element(es["state"].dag(), es["state"]) * qt.ket2dm(es["state"])
231.
232.
233.         else:
234.             dm = qt.ket2dm(states[0][0])
235.             for state in states[1:]:
236.                 if round(state[1], 9) == round(gs_energy, 9):
237.                     dm += qt.ket2dm(state[0])
238.                     #raw_input()
239.
240.         print "Ground state energy: {}".format(gs_energy)
241.         #print """Initial Density matrix:
242.         #        {}""".format(dm)
243.         return dm/dm.tr()
244.
245.     def get_dissipators(self):
246.         dissipators = []
247.         done = []
248.
249.         #here im double counting cause i replace both 0 and 1 with operators on that cluster.
250.         for i in xrange(2**len(self.clusters)):#cycle over all possible ising states
251.             b = format(i, "0{}b".format(len(self.clusters)) )
252.
253.             base = []
254.             for j, s in enumerate(b):
255.                 base += [self.clusters[j].s_dm(int(s)), self.clusters[j].identity_ancilla ]
256.             print b#, base
257.             #now we have a base vector containing this particular b state
258.             for j, c in enumerate(self.clusters):
259.                 q = list(b)
260.                 q[j] = 'd'
261.
262.                 if q not in done:
263.                     done.append(q)
264.                     for n in xrange(c.N + 1):
265.                         m = c.N/2. - n
266.                         Delta_e = c.gamma - 2 * c.beta * m / c.N - 0.5 * sum([self.parameters.J[j,l] * np.sign(0.5 - int(b[l])) for l in xrange(len(self.clusters))])
267.                         vec = list(base)
268.
269.                         vec[2*j], vec[2*j+1] = c.Sp, c.m_dm(n)
270.                         rate = (c.eta * HST(Delta_e))**0.5
271.                         dissipators.append( rate * qt.tensor(vec) )
272.
273.                         vec[2*j], vec[2*j+1] = c.Sm, c.m_dm(n)
274.                         rate = (c.eta * HST(-Delta_e))**0.5
275.                         dissipators.append( rate * qt.Qobj(qt.tensor(vec).full()) )
276.
277.
278.         print len(dissipators)#, dissipators[0]
279.         return dissipators
280.
281.
282.     def prepare_computation(self):
283.         startTime = time.time()
284.         p = self.parameters
285.         self.local_hamiltonians = []
286.         for c in self.clusters:
287.             self.local_hamiltonians.append(c.local_hamiltonian)
288.
289.         h0 = self.get_total_hamiltonian(0)
290.
291.         self.density_matrix = self.get_dm(h0)
292.
293.         self.tlist = np.linspace(0, p.tInterval, num=p.tSteps, endpoint=True)
294.         print "J={} Delta={}".format(p.J, p.Delta)
295.
296.         #print "Total Hamiltonian at t=0:{}".format(h0)
297.
298.     def compute_dynamics(self):
299.
300.         def L(t, args):
301.             L0 = qt.Qobj(qt.liouvillian(self.get_total_hamiltonian(t), []).full())
302.             return L0.data + L1.data
303.
304.         p = self.parameters
305.
306.         dissipators = self.get_dissipators()
307.         print "There are {} dissipators".format(len(dissipators))
308.         L1 = qt.Qobj(qt.liouvillian(self.identity, dissipators).full())
309.
310.         operators_compute = self.get_operators_compute()
311.         try:
312.             #print self.tlist
313.             opts = qt.Odeoptions(nsteps = 10000000, rhs_reuse=True, min_step=0.001, max_step=0, num_cpus = 4, rtol = 0.001, atol = 0.00001, store_states = True)
314.             self.expt_list = qt.mesolve(L, self.density_matrix, self.tlist, [], operators_compute, progress_bar = True, options=opts)
315.
316.         except Exception as e:
317.             print e
318.             self.expt_list = "Error"
319.
320.         if not os.path.exists("archive/annealing_{}".format(self.folder)):
321.             os.makedirs("archive/annealing_{}".format(self.folder))
322.         for i, r in enumerate(self.expt_list.expect):
323.             with open("archive/annealing_{}/expectation_{}".format(self.folder, i), "w") as f:
324.                 for t, res in zip(self.expt_list.times, r):
325.                     f.write("{} {}".format(t, res) + "\n")
326.
327.         with open("archive/annealing_{}/parameters".format(self.folder), "w") as f:
328.             for l in ["{} = {}".format(k, v) for k, v in self.parameters.__class__.__dict__.iteritems() ]:
329.                 f.write(l + "\n")
330.
331.         return self.expt_list
332.
333.
334.
335.
336.
337.
338.     def show_states(self):
339.         self.eigenstates_1 = self.calculate_eigenstates(p.N1,p.b1, p.gamma1, p.alfa1, self.local_hamiltonian_1)
340.         self.eigenstates_2 = self.calculate_eigenstates(p.N2,p.b2, p.gamma2, p.alfa2, self.local_hamiltonian_2)
341.
342.         print "Hamiltonian 1: b= {}  g= {}".format(p.b1, p.gamma1)
343.         print "Eigenstates: (energy)(theorEnergy) State |ps M I>", len(self.eigenstates_1)
344.         for gs in self.eigenstates_1:
345.             print "(%.3f)(%.3f) |%s %+g %g> = %s"%(gs["energy"],gs["theorEnergy"],gs["ps"], gs["M"],gs["I"], gs["representation"])
346.
347.         print "Hamiltonian 2: b= {}  g= {}".format(p.b2, p.gamma2)
348.         print "Eigenstates: (energy)(theorEnergy) State |ps M I>", len(self.eigenstates_2)
349.         for gs in self.eigenstates_2:
350.             print "(%.3f)(%.3f) |%s %+g %g> = %s"%(gs["energy"],gs["theorEnergy"],gs["ps"], gs["M"],gs["I"], gs["representation"])
351.
352.         statesTime = time.time()
353.         print "States and operators got in %f"%(statesTime-startTime)

Paste is for source code and general debugging text.

Login or Register to edit, delete and keep track of your pastes and more.

Recent Pastes

Raw Paste

or to edit or fork this paste. It's free.