2D Hubbard Model
Complete source code can be downloaded from DMRG.
clear;
%%% begin input parameters %%%
t = 1; % nn-hopping amplitude
t_ = -0.2; % nnn-hopping amplitude
U = 8; % onsite repulsion
width = 5; % width of cylinder
height = 4; % height of cylinder
fPBC = true; % periodic boundary conditions along y
sym = 'SU2'; % U(1)_{charge} x SU(2)_{spin}
%sym = 'U1'; % U(1)_{charge} x U(1)_{spin}
D = 2000; % MPS bond dimension
w = ham.get_num_mpo_states( sym, height, fPBC, t_ );
D_prime = floor(D/w);
D_hat = D;
D_tilde = floor(D/10);
num_sweeps = 1000; % number of sweeps
num_davidson = 3; % number of davidson steps
eps_davidson = 1e-10; % threshold for termination of davidson
num_zero = 1e-10; % numerical zero
% pick type of optimization
%opt_enum = 1; % twosite
opt_enum = 2; % shrewd cbe
%%% end input paramters %%%
N = width*height;
if strcmp(sym,'U1')
[F,Z,S,I] = getLocalSpace('FermionS','Acharge,Aspin');
else
[F,Z,S,I] = getLocalSpace('FermionS','Acharge,SU2spin');
end
local_space = {F,Z,S,I};
% initialize tensor network
network = Network;
network.mps = ham.get_prod_state(N,sym);
network.mpo = ham.get_mpo( t, t_, U, sym, width, height, fPBC );
network.model_iden = I.E;
network.D = D;
network.D_prime = D_prime;
network.D_hat = D_hat;
network.D_tilde = D_tilde;
network.num_davidson = num_davidson;
network.eps_davidson = eps_davidson;
network.num_zero = num_zero;
network.num_zero = num_zero;
network.opt_enum = opt_enum;
network.timetracker = TimeTracker;
fprintf("DMRG ground state calculation for the Hubbard model on a cylinder\n\n");
if strcmp(sym,'U1')
fprintf("Acharge,Aspin\n");
else
fprintf("Acharge,SU2spin\n");
end
if fPBC==1
fprintf("periodic boundary conditions in y\n");
else
fprintf("open boundary conditions in y\n");
end
fprintf("t: %.16f\n",t);
fprintf("t_: %.16f\n",t_);
fprintf("U: %.16f\n",U);
fprintf("width: %i\n",width);
fprintf("height: %i\n",height);
fprintf("D: %i\n",D);
fprintf("D_prime: %i\n",D_prime);
fprintf("D_hat: %i\n",D_hat);
fprintf("D_tilde: %i\n",D_tilde);
fprintf("num_sweeps: %i\n",num_sweeps);
fprintf("num_davidson: %.16f\n",num_davidson);
fprintf("eps_davidson: %.16f\n",eps_davidson);
fprintf("num_zero: %i\n",num_zero);
if opt_enum==1
fprintf("twosite opt\n");
end
if opt_enum==2
fprintf("shrewd cbe opt\n");
end
fprintf("\n\n");
% contract from right to left in the beginning
for i=(N:-1:3)
network.env_contract(i,i+1);
end
if opt_enum==1
% sweep and optimize twosite
for i=(1:num_sweeps)
network.twosite_sweep(i);
network.print_time();
pid = sprintf('%i',feature('getpid'));
[tmp mem_usage] = system(['cat /proc/' pid '/status | grep VmRSS']);
fprintf("RSS: %i MB\n", round(str2num(strtrim(extractAfter(extractBefore(mem_usage, ' kB'), ':'))) / 1000));
[tmp mem_usage] = system(['cat /proc/' pid '/status | grep VmHWM']);
fprintf("Peak RSS: %i MB\n\n", round(str2num(strtrim(extractAfter(extractBefore(mem_usage, ' kB'), ':'))) / 1000));
end
else
% sweep and optimize cbe-onesite
for i=(1:num_sweeps)
network.cbe_sweep(i);
network.print_time();
pid = sprintf('%i',feature('getpid'));
[tmp mem_usage] = system(['cat /proc/' pid '/status | grep VmRSS']);
fprintf("RSS: %i MB\n", round(str2num(strtrim(extractAfter(extractBefore(mem_usage, ' kB'), ':'))) / 1000));
[tmp mem_usage] = system(['cat /proc/' pid '/status | grep VmHWM']);
fprintf("Peak RSS: %i MB\n\n", round(str2num(strtrim(extractAfter(extractBefore(mem_usage, ' kB'), ':'))) / 1000));
end
end