gusucode.com > qit_matlab_0.10.0工具箱源码程序 > qit/+markov/@bath/bath.m

    classdef bath < handle
% BATH  Handle class for heat baths.
%
%  Currently only one type of bath is supported, a bosonic
%  canonical ensemble at absolute temperature T, with a
%  single-term coupling to the system.
%
%  The bath spectral density is Ohmic with a cutoff.
%    J(omega) = hbar^2*omega*cut(omega)*heaviside(omega);
%
%  Two types of cutoffs are supported: exponential and sharp.
%    cut(omega) = exp(-omega/omega_c)
%    cut(omega) = heaviside(omega_c - omega)

% gamma(omega) == 2*pi/hbar^2 (J(omega)-J(-omega))(1+n(omega))

% Ville Bergholm 2009-2010


properties (SetAccess = protected)
  % basic parameters
  type   % Bath type. Currently only 'ohmic' is supported.
  omega0 % hbar*omega0 is the unit of energy for all Hamiltonians (omega0 in Hz)
  T      % Absolute temperature of the bath (in K).

  % shorthand
  scale  % dimensionless temperature scaling parameter hbar*omega0/(kB * T)

  % additional parameters with default values
  cut_type  % spectral density cutoff type
  cut_limit % spectral density cutoff angular frequency/omega0
    
  % spectral density
  % J(omega0 * x)/omega0 = \hbar^2 * j(x) * heaviside(x) * cut_func(x);
  j        % spectral density profile
  cut_func % cutoff function

  % lookup tables / function handles for g and s
  % gamma(omega0 * x)/omega0 = g_func(x) * cut_func(x)
  % gamma(omega0 * dh(k))/omega0 = gs_table(1,k)
  %     S(omega0 * dh(k))/omega0 = gs_table(2,k)
  g_func
  g0
  s0
  dH
  gs_table
end  


methods
function b = bath(type, omega0, T)
  % constructor
  %  b = bath(type, omega0, T)
  %
  %  Sets up a descriptor for a heat bath coupled to a quantum system.
    
    
  global qit;

  if (nargin ~= 3)
    error('Usage: bath(type, omega0, T)')
  end

  % basic bath parameters
  b.type   = type;
  b.omega0 = omega0;
  b.T      = T;

  % shorthand
  b.scale = qit.hbar * omega0 / (qit.kB * T);

  switch type
    case 'ohmic'
      % Ohmic bosonic bath, canonical ensemble, single-term coupling

      b.g_func = @(x) 2*pi*x.*(1 +1./(exp(b.scale*x)-1));
      b.g0 = 2*pi/b.scale; % limit of g at x == 0
      b.j = @(x) x;

    otherwise
      error('Unknown bath type.')
  end

  % defaults, can be changed later
  set_cutoff(b, 'sharp', 20);
end


function build_LUT(b)
  % Build a lookup table for the S integral.

  error('unused')

  % TODO justify limits for S lookup
  if (limit > 10)
    temp = logspace(log10(10.2), log10(limit), 50);
    temp = [linspace(0.1, 10, 100), temp]; % sampling is denser near zero, where S changes more rapidly
  else
    temp = linspace(0.1, limit, 10);
  end
  b.dH = [-fliplr(temp), 0, temp]

  b.s_table = [];
  for k=1:length(b.dH)
    b.s_table(k) = S_func(b, b.dH(k));
  end

  b.dH      = [-inf, b.dH, inf];
  b.s_table = [0, b.s_table, 0];

  plot(b.dH, b.s_table, 'k-x');
end


function S = S_func(b, dH)
% Compute S(dH*omega0)/omega0 = P\int_0^\infty dv J(omega0*v)/(hbar^2*omega0) (dH*coth(scale/2 * v) +v)/(dH^2 -v^2)

  ep = 1e-5; % epsilon for Cauchy principal value
  if (abs(dH) <= 1e-8)
    S = b.s0;
  else
    % Cauchy principal value, integrand has simple poles at \nu = \pm dH.
    f = @(nu) b.j(nu) .* b.cut_func(nu) .* (dH*coth(b.scale * nu/2) +nu) ./ (dH^2 -nu.^2);
    S = quad(f, ep, abs(dH)-ep)...
        +quad(f, abs(dH)+ep, 100*b.cut_limit); % FIXME upper limit should be inf, this is arbitrary
  end
end


function set_cutoff(b, type, lim)
  b.cut_type = type;
  b.cut_limit = lim; % == omega_c/omega0

  % update cutoff function (at least Octave uses early binding, so when
  % parameters change we need to redefine it)
  switch b.cut_type
    case 'sharp'
      b.cut_func = @(x) (abs(x) <= b.cut_limit); % Heaviside theta cutoff
    case 'exp'
      b.cut_func = @(x) exp(-abs(x)/b.cut_limit); % exponential cutoff
    otherwise
      error('Unknown cutoff type "%s"', b.cut_type)
  end

  switch b.type
    case 'ohmic'
      b.s0 = -b.cut_limit; % limit of S at dH == 0
  end
  
  % clear lookup tables, since changing the cutoff requires recalc of S
  % start with a single point precomputed
  b.dH = [-inf, 0, inf];
  b.gs_table = [0, b.g0, 0;
		0, b.s0, 0];
end

end

end