PYTHON 6
SpinProva.py Guest on 20th November 2020 05:57:16 PM
  1. # -*- coding: UTF-8 -*-
  2. import sys, time, math
  3. import qutip as qt
  4. import numpy as np
  5. import copy
  6.  
  7. def decompose(N, s=0.5):
  8.         global Ix, Iy, Iz, Isq, Ip
  9.         Ix = sum([nuclearOperator(qt.jmat(s, 'x'), n, N, s) for n in xrange(N)])#sum([nuclearOperator(0.5 * qt.sigmax(), n, N) for n in xrange(N)])
  10.         Iy = sum([nuclearOperator(qt.jmat(s, 'y'), n, N, s) for n in xrange(N)])#sum([nuclearOperator(0.5 * qt.sigmay(), n, N) for n in xrange(N)])
  11.         Iz = sum([nuclearOperator(qt.jmat(s, 'z'), n, N, s) for n in xrange(N)])#sum([nuclearOperator(0.5 * qt.sigmaz(), n, N) for n in xrange(N)])
  12.         Isq = (Ix**2 + Iy**2 + Iz**2).tidyup()
  13.        
  14.         Ip = sum([nuclearOperator(qt.jmat(s, '+'), n, N, s) for n in xrange(N)])#sum([nuclearOperator(qt.sigmap(), n, N) for n in xrange(N)])
  15.         Im = sum([nuclearOperator(qt.jmat(s, '-'), n, N, s) for n in xrange(N)])#sum([nuclearOperator(qt.sigmam(), n, N) for n in xrange(N)])
  16.         ALLSTATES = []
  17.         stuff = qt.simdiag([Iz, Isq])
  18.         tot = int((2*s+1))**N
  19.        
  20.         #print stuff
  21.         p = 0
  22.         for i in xrange(tot):
  23.                 M = round(stuff[0][0][i] + 0,2)#to be sure that 0 is always positive
  24.                 I = round((-1. + (1+ 4* stuff[0][1][i])**0.5)/2. +0 ,2)
  25.                 state = stuff[1][i].tidyup().unit()
  26.                 if M == +I:    
  27.                         p += 1         
  28.                         ALLSTATES.append({'I' : I,'M' : M,'p' : p,'state' : state})
  29.                         #print I,M,p,state
  30.        
  31.         ALLSTATES = orthogonalize(ALLSTATES)
  32.        
  33.         checkOverlap(ALLSTATES)
  34.         #return ALLSTATES
  35.         #print ALLSTATES
  36.         for elem in ALLSTATES:         
  37.                 if elem['M'] - 1 >= -elem['I']:
  38.                         Im2 = qt.Qobj(Im.full())
  39.                         state2 = qt.Qobj(elem['state'].full())
  40.                         elem['state'] = state2
  41.                         newState= (Im2 * state2)
  42.                         newElem = dict(elem)
  43.                         newElem['M'] -= 1
  44.                         #print newState
  45.                         try:   
  46.                                 newElem['state'] = newState.unit().tidyup()
  47.                         except IndexError, ZeroDivisionError:
  48.                                 print elem['state'],elem['I'],elem['M'], newState
  49.                                 pass
  50.                                 #print elem['state'],elem['I'],elem['M'], newState
  51.                                
  52.                         ALLSTATES.append(newElem)
  53.                
  54.  
  55.         ALLSTATES = sorted([el for el in ALLSTATES], key = lambda el: (el['I'], el['M']))
  56.        
  57.         return ALLSTATES
  58. def orthogonalize(allStates):
  59.         def gramSchmidt(states):
  60.                 orthostates = []
  61.                 for v1 in states:
  62.                         for v2 in orthostates:
  63.                                 v1 = v1 - v2.overlap(v1.dag()) * v2
  64.                         orthostates.append(v1.unit().tidyup())
  65.                 return orthostates
  66.  
  67.         orthostates = gramSchmidt([s['state'] for s in allStates])
  68.         for i in xrange(len(allStates)):
  69.                 allStates[i]['state'] = orthostates[i]
  70.         #print allStates
  71.         return allStates
  72.  
  73.  
  74. def represent(allstates):
  75.         tot = len(allstates)
  76.         print 'There are %d nuclear states'%tot
  77.         N = int(math.log(tot,2))
  78.         for elem in allstates:
  79.                 print 'I = %g  M = %g p = %g'%(elem['I'],elem['M'],elem['p'])
  80.                
  81.                 #create a list with all possible permutations of 0 and 1 for each spin
  82.                 spinsList = [[0 for _ in xrange(N)]]
  83.                 for t in xrange(1,tot):
  84.                         binT = bin(t)[2:].zfill(N)
  85.                         spinsList.append([int(binT[i]) for i in xrange(N)])
  86.                 #to each entry in the state correspond a precise spins configuration.
  87.                 for i,d in enumerate(elem['state'].full()):
  88.                         #print elem['state'].full()
  89.                         if d.real != 0:
  90.                                 print '%+.2f %2s'%(d.real , '|'+' '.join([UPARROW if s==0 else DOWNARROW for s in spinsList[i]])+' >'),
  91.                
  92.                 print '\n'
  93.  
  94.  
  95. def nuclearOperator(QObj, n, N, s=1./2):
  96.        
  97.         nucleiList = [qt.qeye(2*s + 1) for _ in xrange(N)]
  98.         nucleiList[n] = QObj
  99.         return qt.tensor(nucleiList)   
  100.  
  101.  
  102. def checkOverlap(allstates):
  103.         for i in xrange(len(allstates)):
  104.                 for j in xrange(i+1,len(allstates)):
  105.                         ol = allstates[i]['state'].overlap(allstates[j]['state']).real
  106.                         if ol < 0.000000001:
  107.                                 ol = 0
  108.                         if ol:
  109.                                 print 'Overlap < %g %g %d|%g %g %d > = %g'%(allstates[i]['I'],allstates[i]['M'],allstates[i]['p'],allstates[j]['I'],allstates[j]['M'],allstates[j]['p'],ol)
  110.                                
  111. UPARROW = u'\u2191'
  112. DOWNARROW = u'\u2193'
  113. CHECKOVERLAP = 1
  114. ORTHOGONAL = 0
  115. if __name__ == '__main__':
  116.         try:   
  117.                 K = int(sys.argv[1])
  118.         except:
  119.                 K = 4
  120.        
  121.         iz = qt.jmat((K/2.), 'z')      
  122.         #print  iz
  123.         for i in xrange(K+1):
  124.                 state = qt.basis(K + 1, i)
  125.                 #print state,iz.matrix_element(state.dag(), state).real
  126.         izc = sum([nuclearOperator(0.5 * qt.sigmaz(), n, K, 0.5) for n in xrange(K)])
  127.  
  128.         for statec in decompose(K):
  129.                 if statec['I'] == K/2.:
  130.                         st = statec['state']
  131.                         #print st,izc.matrix_element(st.dag(), st).real
  132.                
  133.         print [izc.eigenstates()[0][_],

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.