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
  39. def Ifunc(x, tau_lead):
  40.     x = x.real
  41.     threshold = 14
  42.     if x == 0:
  43.         return tau_lead
  44.    
  45.     elif tau_lead == 0:
  46.         if x > 0:
  47.             return 0
  48.         elif x < 0:
  49.             return -x
  50.    
  51.     elif x/tau_lead < -threshold:
  52.         return -x
  53.     elif x/tau_lead > threshold:
  54.         return 0
  55.     return (x / (math.exp(x/tau_lead)-1))
  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.

Raw Paste

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