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

    function [H, dim] = heisenberg(dim, C, J, B)
% HEISENBERG  Heisenberg spin network model
%  [H, dim] = heisenberg(dim, C, J, B)
%
%  Returns the Hamiltonian H for the Heisenberg model, describing a
%  network of n interacting spins in an external magnetic field.
%
%  dim is a dimension vector of the spins, i.e. dim == [2 2 2]
%  would be a system of three spin-1/2's.
%
%  C is the n \times n connection matrix of the spin network, where C(i, j)
%  is the coupling strength between spins i and j. Only the upper triangle is used.
%
%  J defines the form of the spin-spin interaction. It is either a
%  3-vector (homogeneous couplings) or a function J(i, j) returning
%  a 3-vector for site-dependent interactions.
%  Element k of the vector is the coefficient of the Hamiltonian term S_k^{(i)} S_k^{(j)},
%  where S_k^{(i)} is the k-component of the angular momentum of spin i.
%
%  B defines the effective magnetic field the spins locally couple to. It's either
%  a 3-vector (homogeneous field) or a function B(i) that returns a 3-vector for
%  site-dependent field.
%
%  Examples:
%    C = diag(ones(1, n-1), 1)  linear n-spin chain
%    J = [2 2 2]                isotropic Heisenberg coupling
%    J = [2 2 0]                XX+YY coupling
%    J = [0 0 2]                Ising ZZ coupling
%    B = [0 0 1]                homogeneous Z-aligned field
%
%  H =  \sum_{\langle i,j \rangle} \sum_{k = x,y,z} J(i,j)[k] S_k^{(i)} S_k^{(j)}
%      +\sum_i \vec{B}(i) \cdot \vec{S}^{(i)})

% Ville Bergholm 2009-2014


if nargin < 1
    error('Usage: heisenberg(dim, C, J, B)')
end
n = length(dim);

if nargin < 2
    % linear open chain
    C = 'open';
end
if nargin < 3
    % Ising ZZ
    J = [0, 0, 2];
end
if nargin < 4
    % Z field
    B = [0, 0, 1];
end

% simple way of defining some common cases
if ischar(C)
    switch C
      case 'open'
        C = diag(ones(1, n-1), 1);
      case 'closed'
        C = diag(ones(1, n-1), 1);
        C(1, end) = 1;
      otherwise
        error('Unknown topology.')
    end
end

if isa(J, 'function_handle')
    Jfunc = J;
else
    if isvector(J) && length(J) == 3
        Jfunc = @(i, j) C(i, j) * J;
    else
        error('J must be either a 3x1 vector or a function handle.')
    end
end

if isa(B, 'function_handle')
    Bfunc = B;
else
    if isvector(B) && length(B) == 3
        Bfunc = @(i) B;
    else
        error('B must be either a 3x1 vector or a function handle.')
    end
end

% local magnetic field terms
temp = {};
for i = 1:n
    A = angular_momentum(dim(i));  % spin ops
    temp = cat(2, temp, {{cdot(Bfunc(i), A), i}});
end
H = op_list(temp, dim);

% spin-spin couplings: loop over nonzero entries of C
% only use the upper triangle
[a, b] = find(triu(C));
for m = 1:length(a)
    i = a(m);
    j = b(m);
    % spin ops for sites i and j
    Si = angular_momentum(dim(i));
    Sj = angular_momentum(dim(j));
    temp = {};
    % coupling between sites a and b
    c = Jfunc(i, j);
    for k = 1:3
        temp = cat(2, temp, {{c(k) * Si{k}, i; Sj{k}, j}});
        H = H +op_list(temp, dim);
    end
end
end

function res = cdot(v, A)
% Real dot product of a vector with a cell vector of operators.

res = 0;
for k = 1:length(v)
    res = res +v(k) * A{k};
end
end