gusucode.com > qit_matlab_0.10.0工具箱源码程序 > qit/+hamiltonian/spin1d.m

    function varargout = spin1d(dim, J, h, ss)
% SPIN1D  Spin model, 1D spin chain.
%  [H, dim] = spin1d(dim, J, h)               % full Hamiltonian
%  [H, dim] = spin1d(dim, J, h, [first last]) % partial H
%  [S_k, S_{k+1}, h_k, h_{k+1}] = spin1d(dim, J, h, k) % operators
%
%  Returns the Hamiltonian H for an implementation of the
%  Hamiltonian model, describing a one-dimensional chain of spins
%  with nearest-neighbor couplings in an external magnetic
%  field in the z direction.
%
%  The secondary calling syntax returns the coupling (S1, S2) and
%  local (h1, h2) operators just for sites k and k+1 in the chain.
%
%  dim is a vector of the dimensions of the spins, i.e. dim == [2 2 2]
%  would be a chain of three spin-1/2's.
%
%  J defines the coupling strengths in the chain. It's either a
%  upper triangular 3x3 matrix (homogeneous couplings) or a function handle
%  J(site) that returns upper triangular matrices (site-dependent couplings).
%
%  By setting J = [0 0 0;
% 		   0 0 0;
% 		   0 0 a] the Ising model is obtained.
%  Likewise, J =  [a 0 0;
% 		   0 a 0;
% 		   0 0 0] yields the XY model
%
%  h defines the couplings of the spins to the magnetic field and is a 3x1 vector
%  (homogeneous) or a function handle h(site) that returns 
%  3x1 vectors (site-dependent).
%
%  H = \sum_k (\sum{a=x,y,z}\sum{b=x,y,z} J_k(a,b) S_k^a S_k^b + h_k(a)S_k^a)

% Ville Bergholm 2009-2010
% James Whitfield 2010


if (nargin < 3 || nargin > 4)
  error('Usage: spin1d(dim, J, h)')
end

if (isa(J, 'function_handle'))
  Jfunc = J;
else
  if (size(J,1)==size(J,2) && length(J) == 3)
    Jfunc = @(k) J;
  else
    error(['J must be either a 3x3 matrix or a function handle.'])
  end
end

if (isa(h, 'function_handle'))
  hfunc = h;
else
  if (isvector(h))
    hfunc = @(k) h;
  else
    error('h must be either a vector or a function handle.')
  end
end

if (~isvector(dim))
  error('Only 1D spin chains.')
end
n = size(dim);
n = prod(n);


if (nargin == 4)
  if (isscalar(ss))
    % just return the ops for sites ss and ss+1
    S1 = angular_momentum(dim(ss));
    S2 = angular_momentum(dim(ss+1));
    J=Jfunc(ss);
    %matJ=
	%[XX,XY,XZ ]
	%[   YY,YZ ]
	%[      ZZ ]
    h=hfunc(ss);
    %h=
	%[X]
	%[Y]
	%[Z]
    idx=1;
    for comp1=1:3
   	 for comp2=(1+(comp1-1)):3
		if J(comp1,comp2)
			varargout{1}{idx} = J(comp1,comp2) * S1{comp1}; % include J in S1
			varargout{2}{idx} = S2{comp2};
			idx=idx+1;
		end
    	end

	varargout{3}{comp1} = h(comp1) * S1{comp1}; % h1
	varargout{4}{comp1} = h(comp1) * S2{comp1}; % h1
    end

    return
  else
    % return partial H, ss == [start, end] (range of sites)
  end
else
  ss = [1, length(dim)]; % H for entire chain (default)
end


H = sparse(0);
S1 = angular_momentum(dim(ss(1))); % spin ops for first site

rdim = dim(ss(1):ss(2)); % reduced dimension vector
m = 1; % reduced subsystem index to go with rdim

for k=ss(1):ss(2)-1
  % coupling between spins k and k+1: S1 \cdot S2
  S2 = angular_momentum(dim(k+1)); % spin ops
  J=Jfunc(k);
  h=hfunc(k);
  %H = H +op_list({{Jfunc(k,1)*S1{1}, m; S2{1}, m+1}, {Jfunc(k,2)*S1{2}, m; S2{2}, m+1},...
  %                {Jfunc(k,3)*S1{3}, m; S2{3}, m+1}, {hfunc(k)*S1{3}, m}}, rdim);

  H = H+op_list({...
  {J(1,1)*S1{1},m; S2{1},m+1},{J(1,2)*S1{1}, m; S2{2}, m+1},{J(1,3)*S1{1}, m; S2{3}, m+1},...
                              {J(2,2)*S1{2}, m; S2{2}, m+1},{J(2,3)*S1{2}, m; S2{3}, m+1},...
                                                            {J(3,3)*S1{3}, m; S2{3}, m+1},...
  {h(1)*S1{1}, m},            {h(2)*S1{2}, m},              {h(3)*S1{3},m}},rdim);

  % ops for next iteration
  S1 = S2;
  m = m+1;
end

% final term
h=hfunc(ss(2));
H = H +op_list({{h(1) * S1{1}, m},{h(2) * S1{2}, m},{h(3) * S1{3}, m}}, rdim);

varargout{1} = H;
varargout{2} = rdim;