global hadamard cnot cphase sx sy sz rx ry rz N zz xx yy %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %pauli matrices sx=[0 1; 1 0]; sy=[0 -i; i 0]; sz=[1 0; 0 -1]; si=[1 0; 0 1]; pauli={sx,sy,sz}; h=sym('h'); J=sym('J'); K=sym('K'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %two-qubit interaction terms: product space zz=kron(sz,sz); xx=kron(sx,sx); yy=kron(sy,sy); zi=kron(sz,si); iz=kron(si,sz); ii=kron(si,si); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %State Vectors--This will create a column vector of length 16 representing %The 16 possible combinations of up and down psi0=[1; 0]; psi1=[0; 1]; psi00=kron(psi0,psi0); psi01=kron(psi0,psi1); psi10=kron(psi1,psi0); psi11=kron(psi1,psi1); psi0000=kron(psi00,psi00); psi0001=kron(psi00,psi01); psi0010=kron(psi00,psi10); psi0011=kron(psi00,psi11); psi0100=kron(psi01,psi00); psi0101=kron(psi01,psi01); psi0110=kron(psi01,psi10); psi0111=kron(psi01,psi11); psi1000=kron(psi10,psi00); psi1001=kron(psi10,psi01); psi1010=kron(psi10,psi10); psi1011=kron(psi10,psi11); psi1100=kron(psi11,psi00); psi1101=kron(psi11,psi01); psi1110=kron(psi11,psi10); psi1111=kron(psi11,psi11); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Basis State Vectors-This will create the state vectors symbolically %u stands for spin up; d stands for spin down uuuu=sym('uuuu');uuud=sym('uuud');uudu=sym('uudu');uudd=sym('uudd'); uduu=sym('uduu');udud=sym('udud');uddu=sym('uddu');uddd=sym('uddd'); duuu=sym('duuu');duud=sym('duud');dudu=sym('dudu');dudd=sym('dudd'); dduu=sym('dduu');ddud=sym('ddud');dddu=sym('dddu');dddd=sym('dddd'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%Here is the basis of spin up down combinations basis=[uuuu;uuud;uudu;uudd;uduu;udud;uddu;uddd;duuu;duud;dudu;dudd;... dduu;ddud;dddu;dddd]; %Test to see if it works for two spin 1/2 particles %I believe it does s1.s2=kron(xx,ii)+kron(yy,ii)+kron(zz,ii); s2.s3=kron(kron(si,sx),kron(sx,si))+kron(kron(si,sy),kron(sy,si))+kron(kron(si,sz),kron(sz,si)); s3.s4=kron(ii,xx)+kron(ii,yy)+kron(ii,zz); s4.s1=kron(kron(sx,si),kron(si,sx))+kron(kron(sy,si),kron(si,sy))+kron(kron(sz,si),kron(si,sz)); s1.s3=kron(kron(sx,si),kron(sx,si))+kron(kron(sy,si),kron(sy,si))+kron(kron(sz,si),kron(sz,si)); s2.s4=kron(kron(si,sx),kron(si,sx))+kron(kron(si,sy),kron(si,sy))+kron(kron(si,sz),kron(si,sz)); H=5*(s1.s2+s2.s3+s3.s4+s4.s1)+2*(s1.s3+s2.s4); [evec,eval]=eig(H); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Test to see if Hamiltonian com %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Print out eigenvectors and eigenvalues in nice fashion for i=1:16 a(i,1)=eval(i,i); a(i,2)=evec(:,i)'*basis; end disp(' Eigenvalue Eigenvector'); disp(a);