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.

Recent Pastes

Raw Paste

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