/
HF_multi.m
133 lines (107 loc) · 3.71 KB
/
HF_multi.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
% Script for solving the HF equations by energy minimization for the
% four-flavor model with a contact interaction
%% Parameters
clear param;
param.nf = 4; % number of flavors
param.W = 1; % bandwidth
param.U = 1; % interaction strength
param.N = 2000;
param.E0 = 0.7; % DOS parameters
param.Evh = 0.7;
param.Ecut = 1/param.N;
param.nrand = 5; % number of random initial guesses
param.B = 0; % Zeeman field (deactivated)
param.spin = [1,-1,1,-1]';
param.is_lin_dos = 'y'; % if param.is_lin_dos=='y', use analytic expressions for linear DOS
%% Initialize DOS(eps), Ek(eps), n(eps)
eps = -param.W:param.W/param.N:param.W;
if (param.is_lin_dos == 'y') % if not using linear DOS
% define DOS also for case is_lin_dos=='y' for compressibility computation
nu = abs(eps)/param.W^2;
else
% Different DOS models
%nu = abs(eps)/param.W^2;
%nu = 0.1+abs(eps)/param.W^2.*sqrt((abs(eps)-param.Evh).^2+param.Ecut^2).^(-0.25);
%nu = abs(eps)/param.W^2.*sqrt((abs(eps)-param.Evh).^2+param.Ecut^2).^(-0.25);
%nu = abs(eps)/param.W^2;
%nu = 0.1+abs(eps)/param.W^2.*log(sqrt((abs(eps)-param.Evh).^2+param.Ecut^2).^(-1));
%nu = 0.1+abs(eps)/param.W^2;
nu = abs(eps)/param.W^2.*(abs(eps)<=param.E0) ...
+ (param.W-abs(eps))/(param.W-param.E0)*param.E0/param.W^2.*(abs(eps)> param.E0);
end
% Compute n as integral of nu(eps)
param.n = cumtrapz(eps, nu);
param.n = param.n - param.n(param.N+1);
param.n = param.n/param.n(end);
% Normalize nu such that n(eps=W)=1
ntot = trapz(eps, nu);
nu = nu*2/ntot;
% Compute Ek as integral of eps*nu
param.Ek = cumtrapz(eps, eps.*nu);
param.Ek = param.Ek - param.Ek(param.N+1);
[I,J] = meshgrid(1:param.nf, 1:param.nf);
param.utmat = I>J;
%% Loop over mu (chemical potential)
mu = 0:0.02:5;
%U = 0.6:0.005:0.7;
n_min = zeros(param.nf, length(mu));
E_min = zeros(1, length(mu));
lb = [-1,-1,-1,-1]';
ub = [ 1, 1, 1, 1]';
n0 = [1,-1,0,0]';
options = optimoptions('fmincon');
options.Display = 'off';
options.OptimalityTolerance=1e-7;
Etr = zeros(1,param.nrand+1);
ntr = zeros(param.nf,param.nrand+1);
tic
for j = 1:length(mu)
param.mu = mu(j); % chemical potential for this iteration
if (param.is_lin_dos == 'y')
f = @(x)En1n2_lin(x,param);
else
f = @(x)En1n2(x,param); % Define f = En1n2 as anonymous function
end
[ntr(:,1),Etr(1)] = fmincon(f,n0,[],[],[],[],lb,ub,[],options);
% Random initial conditions
n0tr = 2*(rand(param.nf,param.nrand)-0.5);
for jr=1:(param.nrand+1)
if (jr==1)
n_init = n0;
else
n_init = n0tr(:,jr-1);
end
[ntr(:,jr),Etr(jr)] = fmincon(f,n_init,[],[],[],[],lb,ub,[],options);
end
[~,I] = min(Etr);
n_min(:,j) = ntr(:,I);
E_min(j) = Etr(I);
n0 = n_min(:,j);
if (mod(j,50)==0)
fprintf('mu = %.4f\n',mu(j))
end
end
n_min = sort(n_min,1);
n_tot = sum(n_min,1);
%% Compute compressibility
nu_a = interp1(param.n,nu,n_min(:)).*(abs(n_min(:))<0.999);
nu_a = reshape(nu_a, size(n_min,1), size(n_min,2));
nu_bar = sum(nu_a./(1-param.U*nu_a),1);
dndm = nu_bar./(1 + param.U*nu_bar);
%% Plot the results
figure
subplot(2,2,1);
plot(eps,nu,'linewidth',1);
prettyfig; xlabel('$$\varepsilon/W$$'); ylabel('$$\nu$$');
subplot(2,2,2);
plot(mu, n_min,'linewidth',1)
prettyfig; xlabel('$$\mu$$'); ylabel('$$n_\alpha$$');
subplot(2,2,3);
plot(n_tot,mu,'linewidth',1);
prettyfig; xlabel('$$n$$'); ylabel('$$\mu$$');
subplot(2,2,4);
[n_flat, dmdn_flat] = flatten_dmdn(n_tot, mu, dndm); % Corrected dmdn for first order transitions
plot(n_flat, dmdn_flat,'linewidth',1)
prettyfig; xlabel('$$n$$'); ylabel('$$d \mu/dn$$');
axis([0 4 0 3])
toc