Дипломная работа: Численный расчёт фазовой диаграммы двух связанных цепочек холодных атомов в рамках модели Хаббарда с притяжением

Внимание! Если размещение файла нарушает Ваши авторские права, то обязательно сообщите нам

В работе методами численного моделирования проверен ряд теоретических результатов (выражение для энергии связи димера и уравнения для некоторых фазовых границ).

В работе получено, что топология фазовых диаграмм двух и более связанных цепочек атомов качественно не меняется при увеличении числа цепочек, но качественно отличается от топологии диаграммы одной цепочки. Этот результат получен в настоящей работе впервые и представляет реальную научную ценность.

Стоит отметить, что в ходе выполнения работы было разработано программное обеспечение, которое будет использовано в дальнейшем при обучении студентов магистерской програмы “Материалы, приборы, нанотехнологии”.

Библиографический список

1. Кашурников В.А., Красавин А.В., Вычислительные методы в квантовой физике: Учебное пособие. М.: МИФИ, 200. - 412с.

2. Hubbard J., Electron correlations in narrow energy bands-3, 4, "Proc. Roy. Soc. A", 1963, v. 276, p. 238; 1964, v. 281, p. 401; 1965, v. 285, p. 541

3. F. Heidrich-Meisner, G. Orso, and A.E. Feiguin “Phase separation of trapped spin-imbalanced Fermi gases in one-dimensional optical lattices”, Phys. Rev. A 81, 053602

4. A. S. Holevo, Statistical definition of observable and the structure of statistical models. Rep. Math. Phys., 1985, v.22, N3.

5. Ландау Л.Д., Лифшиц Е.М., Краткий курс теоретической физики: Квантовая физика. М., 1972, .-368с.

Приложение 1

Текст программы расчета энергии спиновом кластера

% Main module

t_init = cputime;

EnterData %preparing data for the solver

NewVariables %declaring new variables

Basis %generating basis

m_Ham = Hamiltonian;

eig(m_Ham)

t_fin = cputime;

fprintf(1,'Calculation time, min: %.0f\n',(t_fin-t_init)/60)

% Result of a-operator act on a state vector

function f = a(type,number,state_vect)

%Arguments:_________________________________________________________

%type - type of operator: 1 - creation,-1 - annihilation;

%number - number of node state (to create or annihilate);

%state_vect - state vector (like [10100101110]).

%Returning value: new state vector.

%___________________________________________________________________

N = size(state_vect,2); %size of the state vector

if iscorr(state_vect)

if appzer(type,number,state_vect) %apparent zeros

f = zeros(1,N);

else %the main procedure

f = (-1)^numofones(number,state_vect)*changestate(type,number,state_vect);

end

else

f = zeros(1,N);

disp('Incorrect format of the state vector!')

end

% Checking format of the state vector

function f = iscorr(state_vect)

%Arguments:_________________________________________________________

%state_vect - state vector (like [10100101110]).

%Returning value: 1 - if the state vector is correct (contains just ones and zeros), 0 - otherwise.

%___________________________________________________________________

N = size(state_vect,2); %size of the state vector

f = 1;

for i=1:N

if state_vect(i)~=0 && state_vect(i)~=1 && state_vect(i)~=-1

f = 0;

break

end

end

% Apparent zero

function f = appzer(type,number,state_vect)

%Arguments:_________________________________________________________

%type - type of operator: 1 - creation,-1 - annihilation;

%number - number of node state (to create or annihilate);

%state_vect - state vector (like [10100101110]).

%Returning value: 1 - if resulting vector is zero, 0 - otherwise.

%___________________________________________________________________

f = 0;

switch type

case 1

if state_vect(number) == 1

f = 1;

end

case -1

if state_vect(number) == 0

f = 1;

end

otherwise

disp('Incorrect operator type!')

f = 1;

end

% Changing one state

function f = changestate(type,number,state_vect)

%Arguments__________________________________________________________

%type - type of operator: 1 - creation,-1 - annihilation;

%number - number of node state (to create or annihilate);

%state_vect - state vector.

%Returning value: state_vect with 1 state changed.

%___________________________________________________________________

f = state_vect;

switch type

case 1 %create an electron

if sov(state_vect) >= 0

f(number) = 1;

else f(number) = -1;

end

case -1 %annihilate an electron

f(number) = 0;

otherwise

disp('Incorrect operator type!')

f = 0*state_vect;

end

% Number of ones before given number

function f = numofones(number,state_vect)

%Arguments:_________________________________________________________

%number - number of node state;

%state_vect - state vector.

%Returning value: number of ones before given number.

%___________________________________________________________________

if number == 1

f = 0;

else

f = sum(state_vect(1:number-1));

end

% Result of a_plus_i_a_j-operator act on a state vector

function v_f = a_plus_i_a_j(i,j,v_F_q);

%Arguments:________________________________________________________________________________

%i - number of a_plus-operator;

%j - number of a_minus-operator;

%v_F_q - state vector (without spin, only positive components).

%Returning value: new state vector.

m = min(i,j); n = max(i,j);

k_perm = (-1)^sum(v_F_q(m:n-1));

if i < j

k_sign = 1;

else

k_sign = -1;

end

%a_j|v_F_q>:

k_zero = 1;

if v_F_q(j)

v_F_q(j) = 0;

else

k_zero = 0;

end

%a_plus_i|v_F_q>:

if v_F_q(i)

k_zero = 0;

else

v_F_q(i) = 1;

end

v_f = k_perm*k_zero*k_sign*v_F_q;

% Generating basis

function Basis

global m_basis N_nodes N_part_up N_part_dwn

m_basis_up = GenBasis(N_nodes,N_part_up); %basis functions for spin-up electrons

N_up = size(m_basis_up,1); %number of them

m_basis_dwn = GenBasis(N_nodes,N_part_dwn); %... spin-down electrons

N_dwn = size(m_basis_dwn,1); %number of them

%General basis:

N_set = N_up*N_dwn;

k = 1;

m_basis = zeros(N_set,2*N_nodes);

for i=1:N_up

for j=1:N_dwn

m_basis(k,:) = [m_basis_up(i,:),m_basis_dwn(j,:)];

k = k+1;

end

end

%Entering all parameters from outer files

%Model parameters:____________________________________________________

global eps0 t fl_bc N_nodes N_part_up N_part_dwn

load 'ModParams.dat' -ascii

eps0 = ModParams(1); t = ModParams(2); fl_bc = ModParams(3);

N_nodes = ModParams(4); N_part_up = ModParams(5); N_part_dwn = ModParams(6);

clear ModParams

% Generating the basis set

function m_f = GenBasis(N_nodes,N_part)

%Arguments:_________________________________________________________

%N_nodes - number of nodes;

%N_part - number of particles.

%Returning value: m_f - matrix of basis set.

%__________________________________________________________________

N_basis = nchoosek(N_nodes,N_part);

m_basis1 = nchoosek(1:N_nodes,N_part);

m_f = zeros(N_basis,N_nodes);

for i=1:N_basis

for j=1:N_part

m_f(i,m_basis1(i,j)) = 1;

end

end

% Hamiltonian matrix

function m_f = Hamiltonian

%Arguments:_________________________________________________________

% none

%Returning value: Hamiltonian matrix.

%___________________________________________________________________

global m_basis N_nodes

N = size(m_basis,1); %dimension of state space

m_f = Ham_pot(N) + Ham_kin(N,m_basis(:,1:N_nodes),m_basis(:,N_nodes+1:2*N_nodes)) + ...

Ham_kin(N,m_basis(:,N_nodes+1:2*N_nodes),m_basis(:,1:N_nodes));

% Hamiltonian matrix (potential part)

function m_f = Ham_pot(N)

%Arguments:________________________________________________________________________________

% N - size of the matrix

%Returning value: potential part of Hamiltonian matrix.

%__________________________________________________________________________________________

global eps0 t m_basis N_nodes

m_f = zeros(N);

for i=1:N

state_vect_up = m_basis(i,1:N_nodes); state_vect_dwn = m_basis(i,N_nodes+1:2*N_nodes);

m_f(i,i) = eps0*state_vect_up*state_vect_dwn';

end

% Hamiltonian matrix (kinetic part)

function m_f = Ham_kin(N,m_basis,m_basis_aux)

%Arguments:_________________________________________________________

%N - size of the matrix;

%m_basis - basis matrix (spin up), m_basis_aux - basis matrix (spin down) or vice-versa.

%Returning value: kinetic part of Hamiltonian matrix.

%__________________________________________________________________________________________

global eps0 t fl_bc

m_f = zeros(N);

N_nodes = size(m_basis,2); %number of nodes

arr_trr1(1:N_nodes) = 0; arr_trr2(1:N_nodes) = 0; %arrays for scalar products

for p=1:N

for q=1:N

if p~=q %just non-diagonal elements

v_F_p = m_basis(p,:); v_F_q = m_basis(q,:);

for i=1:N_nodes %sum over the nodes_____________________________________________________

%A+(i)*A-(i-1):___________________________________________________

%Determining number of A-(i-1) operator:

j = i - 1;

fl_calc = 1; %flag indicator: to calc or not to calc the matrix element

if j == 0

switch fl_bc

case 0 %the chain is not closed

fl_calc = 0;

arr_trr1(i) = 0;

case 1 %the chaine is closed

j = N_nodes;

otherwise

disp('Wrong type of boundary conditions!')

end

end

%Calculating the first item of the matrix element:

if fl_calc

v_tr = a_plus_i_a_j(i,j,v_F_q);

arr_trr1(i) = scpr([v_F_p,m_basis_aux(p,:)],[v_tr,m_basis_aux(q,:)]);

end

%A+(i)*A-(i+1):_______________________________________________________________

%Determining number of A-(i-1) operator:

j = i + 1;

fl_calc = 1;

if j == N_nodes + 1

switch fl_bc

case 0 %the chaine is not closed

fl_calc = 0;

arr_trr2(i) = 0;

case 1 %the chaine is closed

j = 1;

otherwise

disp('Wrong type of boundary conditions!')

end

end

%Calculating the second item of the matrix element:

if fl_calc

v_tr = a_plus_i_a_j(i,j,v_F_q);

arr_trr2(i) = scpr([v_F_p,m_basis_aux(p,:)],[v_tr,m_basis_aux(q,:)]);

end

%...A+(i)*A-(i+1)_______________________________________________________________________

end %sum over the nodes_________________________________________________________________

m_f(p,q) = t*( sum(arr_trr1) + sum(arr_trr2) );

end

end

end

load 'H_example.dat' -ascii

m_Hex = H_example;

clear H_example

% Main module

t_init = cputime;

EnterData %preparing data for the solver

NewVariables %declaring new variables

Basis %generating basis

m_Ham = Hamiltonian;

eig(m_Ham)

t_fin = cputime;

fprintf(1,'Calculation time, min: %.0f\n',(t_fin-t_init)/60)

% Result of n-operator act on a state vector

function f = n(state_vect)

%Arguments:________________________________________________________________________________

%state_vect - state vector.

%Returning value: number of particles in this state.

%__________________________________________________________________________________________

f = sum(abs(state_vect));

% Eigenvalue of n_spin_down-operator on a state vector

function f = n_spin_dwn_i(state_vect_spin,i)

%Arguments:________________________________________________________________________________

%state_vect_spin - state vector with spin;

%i - number of node.

%Returning value: number of particles in spin up mode in this state.

%__________________________________________________________________________________________

global N_nodes

state_vect_dwn = state_vect_spin(N_nodes+1:2*N_nodes);

f = abs(state_vect_up(i));

% Eigenvalue of n_spin_up-operator act on a state vector

function f = n_spin_up_i(state_vect_spin,i)

%Arguments:________________________________________________________________________________

%state_vect_spin - state vector with spin;

%i - number of node.

%Returning value: number of particles in spin up mode in this state.

%___________________________________________________________________

global N_nodes

state_vect_up = state_vect_spin(1:N_nodes);

f = abs(state_vect_up(i));

%New parameters:

% m_basis - basis of wave functions

global m_basis

% Scalar product of two state vectors

function f = scpr(state_vect1,state_vect2)

%Arguments:________________________________________________________________________________

%state_vect1,2 - state vectors.

%Returning value: scalar product.

%__________________________________________________________________________________________

N1 = size(state_vect1,2); %size of the first state vector

N2 = size(state_vect2,2); %size of the second state vector

if N1 ~= N2

disp('Different sizes of state vectors!')

f = 0;

else %main procedure

v1 = abs(state_vect1); v2 = abs(state_vect2);

if sum(v1) == 0 && sum(v2) == 0 %if one of them =0

f = 0

else %procedure for nonzero vectors

f = 1;

for i = 1:N1

if v1(i) ~= v2(i)

f = 0;

break

end

end

if f

f = f*sov(state_vect1)*sov(state_vect2);

end

end

end

% Sign of vector

function f = sov(state_vect)

%Arguments:________________________________________________________________________________

%state_vect - state vector.

%Returning value: sign of first nonzero component of the vector.

%__________________________________________________________________________________________

N = size(state_vect,2); %size of the first state vector

f = 1;

for i = 1:N

if state_vect(i) ~= 0

f = sign(state_vect(i));

break

end

end