%-------------------------------------------------------------------------------- % Compute energies and for N=(1:30) and compare with ab initio % results of: % - Chang and Bertsch, PRA 76, 021603(R) (2007) [arXiv:physics/0703190v1] % - von Stecher, Greene and Blume, PRL 99, 233201 (2007) [arXiv:0708.2734v2] % and arXiv:0801.2747v1 % % % SLDA with Bessel-DVR method, Littlejohn and Cargo, J.Chem.Phys. 117, 27 (2002) % Broyden updates on all sites at once % fixed particle numbers for both spin up and down species % smoothing around Ec and finite temperature % % Units: m=hbar=omega=1 since a gas at unitarity in a harmonic trap has % universal properties as well and m, hbar and omega factor out of the % problem, see Bulgac, PRA, 040502(R) (2007) % % This is the acutal program used for calculations in the above publication % % %-------------------------------------------------------------------------------- clear format short g amix = 1.0; % weight in the usual iterative procedure with linear mixing bmix = 1.0; % size of dx step with Broyden cmix = 1.0; % initial Jacobian convergence_limit = 1.0e-9; deps = 1.0e-15; % essentially defines machine precisions ddmu = 1.0e-6; part_c = 0; ener_c = 0; nlpr = 101; nplot = 101; Nmax = 100; Niter = 0*Nmax; Kmax = sqrt(2*Nmax+3)+3; % Rmax = Kmax; Rmax = 9.0; %--------------------------------------------------------------------------------- % The energy cutoff Ec is made smooth with a Fermi function % % A small finite temperature is introduced in order to faciltitate convergence % especially in the case of odd or polarized systems. The actual value of % T is so small as to be insignificant physically and it ca n be put T=0 % at the end of calculations %--------------------------------------------------------------------------------- Ec = 0.95*(Nmax+1.5); dEc = Ec/50; Temp = 0.025; four_pi = 4*pi; %------------------------------------------------ % SLDA parameters %------------------------------------------------ % %------------------------------------------------ % % the original Carlson et al estimates % %------------------------------------------------ % xi = 0.44; % zeta = xi; % eta = 0.486; % alpha = 1.1175; % beta = -0.51955; % igamma = -0.095509; %--------------------------------------------------- % these are likely the best estimates %--------------------------------------------------- xi = 0.42; zeta = xi; eta = 0.504; alpha = 1.14; beta = -0.55269; igamma = -0.090585; bbar = beta - eta^2*(3*pi^2)^(2/3)*igamma/6; %--------------------------------------------------- %------------------------------------------------------------- % here I construct various things needed for the DVR method %------------------------------------------------------------- % determine the zeros of the Bessel functions j0(z) and j1(z) %------------------------------------------------------------- Nzeros = 500; nn = (1:Nzeros); Z0 = pi*nn; Z1 = pi*(nn + 0.5); for in = 1:20 Z1 = nn*pi + atan(Z1); end %------------------------------------------------------------- iz0 = find( Z0 <= Rmax*Kmax); iz1 = find( Z1 <= Rmax*Kmax); %---------------------------------------------------------------------- % Here I construct kinetic energy matrix for l=0 and l=1, % see Littlejohn and Cargo, J.Chem.Phys. 117, 27 (2002) for Bessel-DVR. % I am using the l=0 lattice points for even ls and l=1 for odd ones, % see Nygaard, Bruun, Schneider, Clark and Feder, PRA 69, 053622 (2004) % [arXiv:cond-mat/0312258v1]. I have also checked rather extensively % on a number of exact solutions. %---------------------------------------------------------------------- for l0 = 0:1 if l0 == 0 n0 = iz0'; z0 = Z0(n0)'; % zeros of spherical Bessel function nu0 = 1/2; [m1,m2] = meshgrid(n0,n0); [zz1, zz2] = meshgrid(z0,z0); t0 = (1+2*(nu0^2-1)./z0.^2)/3; t1 = 8*(-1).^(m1-m2).*zz1.*zz2./( (zz1.^2-zz2.^2).^2 + eps ); t1 = t1-diag(diag(t1)); T0 = Kmax^2*(diag(t0) + t1)/2*alpha; zmax0 = z0(end); r0 = z0/Kmax; r0_2 = r0.^2; mm0 = size(n0,1); elseif l0 == 1 n1 = iz1'; z1 = Z1(n1)'; nu0 = 3/2; [m1,m2] = meshgrid(n1,n1); [zz1, zz2] = meshgrid(z1,z1); t0 = (1+2*(nu0^2-1)./z1.^2)/3; t1 = 8*(-1).^(m1-m2).*zz1.*zz2./( (zz1.^2-zz2.^2).^2 + eps ); t1 = t1-diag(diag(t1)); T1 = Kmax^2*(diag(t0) + t1)/2*alpha; zmax1 = z1(end); r1 = z1/Kmax; r1_2 = r1.^2; mm1 = size(n1,1); end end %-------------------------------------------------------------------- %================================================================= % convert from one coordinate set f(r_a) to another f(r_b) % calculate weights %================================================================= U1_0 = zeros(mm1,mm0); U0_1 = zeros(mm0,mm1); C0 = zeros(mm0,1); % 1/F0(z0) C1 = zeros(mm1,1); % 1/F1(z1) for i0 = 1:mm0 a = cos(z0(i0))/sqrt(z0(i0)); b = sin(z1)./sqrt(z1); U1_0(:,i0) = 2*sqrt(z0(i0)*z1)./(z1.^2-z0(i0).^2).*b/a; C0(i0) = sqrt(pi/Kmax); end C0_2 = C0.^2; for i1 = 1:mm1 a = sin(z1(i1))/sqrt(z1(i1)); b = -cos(z0)./sqrt(z0); U0_1(:,i1) = 2*sqrt(z1(i1)*z0)./(z0.^2-z1(i1).^2).*b/a; C1(i1) = sqrt(pi/Kmax/(sin(z1(i1)))^2); end C1_2 = C1.^2; %----------------------------------------------------------------------- % clear Nzeros nn Z0 Z1 m1 m2 zz1 zz2 iz0 iz1 t0 t1 nu0 a b %----------------------------------------------------------------------- % this concludes the determination of various quantities needed for DVR %----------------------------------------------------------------------- Number_Density = []; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for N = 2:30 % total particle number tic N_a = ceil(N/2); % number of spin-up particles N_b = floor(N/2); % number of spin-down particles %---------------------------------------------------------------------- % Start iterations % % The intial guesses for densities, fields and chemical potentials % are rather crude. Note that the particle numbers obtained from these % initial guess for number densities are rather far from the desired % values (almost by a factor of two sometimes), but that has little % effect on the convergence % % %---------------------------------------------------------------------- mu = (3*N)^(1/3)*sqrt(xi); mu_a = mu; % chemical potential for spin-up mu_b = mu; % chemical potential for spin-down mu0 = mu; E0 = (3*N)^(4/3)/4*sqrt(xi); ir_0 = ( (2*mu - r0_2)>0 ); rho_0_a = ( (2*mu - r0_2)/(alpha*(1+bbar)) ).^1.5/(6*pi^2).*ir_0; rho_0_b = ( (2*mu - r0_2)/(alpha*(1+bbar)) ).^1.5/(6*pi^2).*ir_0; rho_0 = rho_0_a + rho_0_b; D_0 = eta*(3*pi^2*rho_0).^(2/3)/2.*ir_0; V_0 = bbar*(3*pi^2*rho_0).^(2/3)/2.*ir_0; ir_1 = ( (2*mu - r1_2)>0 ); rho_1_a = ( (2*mu - r1_2)/(alpha*(1+bbar)) ).^1.5/(6*pi^2).*ir_1; rho_1_b = ( (2*mu - r1_2)/(alpha*(1+bbar)) ).^1.5/(6*pi^2).*ir_1; rho_1 = rho_1_a + rho_1_b; D_1 = eta*(3*pi^2*rho_1).^(2/3)/2.*ir_1; V_1 = bbar*(3*pi^2*rho_1).^(2/3)/2.*ir_1; x0 = [V_0;D_0;V_1;D_1;mu_a;mu_b]; x1 = x0; dx = x0; G0 = x0; G1 = x0; dG = x0; K0 = cmix*eye(size(x0,1)); K1 = K0; convergence = 1; iter = 0; while convergence > convergence_limit %-------------------------------------- iter = iter + 1; Etot = 0; Eqs = []; Ll = []; % nocc_a = []; % nocc_b = []; rho_0_a = zeros(mm0,1); rho_0_b = zeros(mm0,1); kappa_0 = zeros(mm0,1); rho_1_a = zeros(mm1,1); rho_1_b = zeros(mm1,1); kappa_1 = zeros(mm1,1); %-------------------------------------- for l=0:Nmax % angular momentum l0 = mod(l,2); if l0 == 0 ll = alpha*(l*(l+1) - l0*(l0+1))/2; Vl = ll./r0_2; Hh_a = T0 + diag(Vl + V_0 + r0_2/2-mu_b); Hh_b = T0 + diag(Vl + V_0 + r0_2/2-mu_a); Dh = diag(D_0); HH = [[Hh_a,Dh];[Dh,-Hh_b]]; elseif l0 == 1 ll = alpha*(l*(l+1) - l0*(l0+1))/2; Vl = ll./r1_2; Hh_a = T1 + diag(Vl + V_1 + r1_2/2-mu_b); Hh_b = T1 + diag(Vl + V_1 + r1_2/2-mu_a); Dh = diag(D_1); HH = [[Hh_a,Dh];[Dh,-Hh_b]]; end [phi, eigen] = eig(HH); ener = diag(eigen); al = (2*l+1)/four_pi; for ie = 1:size(ener,1) en_t = ener(ie)/(Temp+eps); if en_t > 40 fe = 0; elseif en_t < -40 fe = 1; else fe = 1/(1+exp(en_t)); end en_c = (abs(ener(ie))-Ec)/dEc; if en_c > 40 fc = 0; elseif en_c < -40 fc = 1; else fc = 1/(1+exp(en_c)); end if fc > 0 if l0 == 0 ur0 = phi( 1: mm0,ie)./C0; vr0 = phi(mm0+1:2*mm0,ie)./C0; vr1 = U1_0*vr0; ur1 = U1_0*ur0; elseif l0 == 1 ur1 = phi( 1: mm1,ie)./C1; vr1 = phi(mm1+1:2*mm1,ie)./C1; vr0 = U0_1*vr1; ur0 = U0_1*ur1; end rho_0_a = rho_0_a + (1-fe)*fc*al*vr0.*vr0; rho_0_b = rho_0_b + fe*fc*al*ur0.*ur0; kappa_0 = kappa_0 + (1-2*fe)*fc*al*vr0.*ur0; rho_1_a = rho_1_a + (1-fe)*fc*al*vr1.*vr1; rho_1_b = rho_1_b + fe*fc*al*ur1.*ur1; kappa_1 = kappa_1 + (1-2*fe)*fc*al*vr1.*ur1; ev = vr0.^2.*(mu_a-ener(ie)-V_0) + vr0.*ur0.*D_0; eu = ur0.^2.*(mu_b+ener(ie)-V_0) - vr0.*ur0.*D_0; Etot = Etot + four_pi*al*sum(((1-fe)*ev+fe*eu).*C0_2)*fc; end % nocc_a = [nocc_a, (1-fe)*fc]; % nocc_b = [nocc_b, fe*fc]; end Eqs = [Eqs, ener']; Ll = [Ll, l*ones(size(ener))']; end % ------------ end l loop rho_0_a = rho_0_a./r0_2; rho_0_b = rho_0_b./r0_2; rho_0 = rho_0_a + rho_0_b; kappa_0 = kappa_0./r0_2/2; rho_1_a = rho_1_a./r1_2; rho_1_b = rho_1_b./r1_2; rho_1 = rho_1_a + rho_1_b; kappa_1 = kappa_1./r1_2/2; %-------------------------------------------------------------- % l=0 sites %-------------------------------------------------------------- k0 = sqrt(2*(mu - r0_2/2 - V_0)/alpha); kc = sqrt(2*(Ec + mu - r0_2/2 - V_0)/alpha); Lc = real( (kc-k0/2.*log((kc+k0)./(kc-k0)))/(2*pi^2*alpha) ); g_eff = 1./( rho_0.^(1/3)*igamma - Lc + eps); last_corr = 1 - D_0.^2/(6*pi^2*alpha^2) ... .*real( log((kc+k0)./(kc-k0))./(k0.*rho_0+eps) ); D_0 = -g_eff.*kappa_0; % pairing field V_0 = beta*(3*pi^2*rho_0).^(2/3)/2 - ... D_0.^2*igamma./(3*(rho_0+eps).^(2/3)); % meanfield V_0 = V_0./last_corr; %------------------------------------------------------------------------ Etot = Etot ... + 0.3*beta*(3*pi^2)^(2/3)*four_pi*sum(r0_2.*rho_0.^(5/3).*C0_2) ... - four_pi*sum(r0_2.*D_0.*kappa_0.*C0_2); %------------------------------------------------------------------------ %-------------------------------------------------------------- % l=1 sites %-------------------------------------------------------------- k0 = sqrt(2*(mu - r1_2/2 - V_1)/alpha); kc = sqrt(2*(Ec + mu - r1_2/2 - V_1)/alpha); Lc = real( (kc-k0/2.*log((kc+k0)./(kc-k0)))/(2*pi^2*alpha) ); g_eff = 1./( rho_1.^(1/3)*igamma - Lc + eps); last_corr = 1 - D_1.^2/(6*pi^2*alpha^2) ... .*real( log((kc+k0)./(kc-k0))./(k0.*rho_1+eps) ); D_1 = -g_eff.*kappa_1; % pairing field V_1 = beta*(3*pi^2*rho_1).^(2/3)/2 - ... D_1.^2*igamma./(3*(rho_1+eps).^(2/3)); % meanfield V_1 = V_1./last_corr; %------------------------------------------------------------------- N0_a = four_pi*sum(r0.^2.*C0.^2.*rho_0_a); N0_b = four_pi*sum(r0.^2.*C0.^2.*rho_0_b); R2_0 = four_pi*sum(r0.^4.*C0.^2.*rho_0); %------------------------------------------------------------------- % This is where Broyden method is applied to perform the update % Notice that the meanfield (V_0,V_1), the pairing field (D_0, D_1) % and the chemical potentials (mu_a, mu_b)are updated, until % convergence is acheived for the desired values of specific particle % numbers (N_a, N_b) and accuracy. The first update (iter==1) % is simple, apart from weight, which is of little relevance. % Notice also that one can run the program with the simple linear % mixing by choosing Niter >0 accordingly. %------------------------------------------------------------------- if iter == 1 G0 = x0 - [V_0;D_0;V_1;D_1; ... mu_a+mu_a*(N_a/N0_a-1);mu_b+mu_b*(N_b/N0_b-1)]; % G0 = x0 - [V_0;D_0;V_1;D_1; ... % mu_a+mu0*(N_a/N0_a-1);mu_b+mu0*(N_b/N0_b-1)]; x1 = x0 - amix*G0; dx = x1 - x0; x0 = x1; elseif iter > 1 G1 = x0 - [V_0;D_0;V_1;D_1; ... mu_a+mu_a*(N_a/N0_a-1);mu_b+mu_b*(N_b/N0_b-1)]; % G1 = x0 - [V_0;D_0;V_1;D_1; ... % mu_a+mu0*(N_a/N0_a-1);mu_b+mu0*(N_b/N0_b-1)]; dG = G1 - G0; ket = dx - K0*dG; ket = ket.*(abs(ket)>deps*(abs(G0)+abs(G1))); bra = dx'*K0; inorm = bmix/(bra*dG); K1 = K0 + ket*bra*inorm; if iter < Niter x1 = x0 - amix* G1; % linear mixing else x1 = x0 - K1*G1; end K0 = K1; dx = x1 - x0; x0 = x1; G0 = G1; end V_0 = x0( 1: mm0,1); D_0 = x0( mm0 +1:2*mm0,1); V_1 = x0(2*mm0 +1:2*mm0+ mm1,1); D_1 = x0(2*mm0+mm1+1:2*mm0+2*mm1,1); mu_a = x0(end-1); mu_b = x0(end ); mu = (mu_a+mu_b)/2; dmu = (mu_a-mu_b)/2; convergence = max(abs([G0;dx])); end % --------- end iterations [N, toc, iter, Nmax, Etot, R2_0, mu, dmu*(abs(dmu)>ddmu)] En(N) = Etot; R2(N) = R2_0; Nn(N) = N; Number_Density = [Number_Density, N*ones(size(r0)), r0, rho_0_a, rho_0_b]; end % loop on particle number %------------------------------------------------------------------------ % The rest of the program is for plotting and comparing SLDA results % ab initio results. %------------------------------------------------------------------------ En(1) = 1.37; R2(1) = 1.37; Nn(1) = 1; %-------------------------------------------------------------------------- % Chang and Bertsch, PRA 76, 021603(R) (2007) [arXiv:physics/0703190v1] %-------------------------------------------------------------------------- np_0 = (1:22); en_0 = [1.5, 2.01, 4.28, 5.1, 7.6, 8.7, 11.3, 12.6, 15.6, 17.2, 19.9, 21.5, ... 25.2, 26.6, 30.0, 31.9, 35.4, 37.4, 41.1, 43.2, 46.9, 49.3]; den_0 = [0, 0.02, 0.04, 0.1, 0.1, 0.1, 0.1, 0.1, 0.2, 0.2, 0.2, 0.3, 0.2, ... 0.4, 0.1, 0.3, 0.2, 0.3, 0.3, 0.4, 0.2, 0.1]; np_v = [2, 3, 8, 14, 20]; en_v = [ 1.9463, 4.1944, 11.914, 25.965, 40.765]; en_vv= (3*np_v).^(4/3)/4*sqrt(xi); %-------------------------------------------------------------------------- in = 0; for N = 3:2:21 in = in + 1; Gap_slda(in) = En(N) -(En(N+1) + En(N-1))/2; Gap_gfmc(in) = en_0(N)-(en_0(N+1) + en_0(N-1))/2; end %-------------------------------------------------------------------------- % J. von Stecher, C.H. Greene, and D. Blume, % PRL 99, 233201 (2007) [arXiv:0708.2734v2] and arXiv:0801.2747v1 % NB There might be slight differences between the two works in the % last digits. %-------------------------------------------------------------------------- doerte = [ 1 1.5 0.0; 2 2.002 0.0; 3 4.281 0.004; 4 5.051 0.009; 5 7.610 0.01; 6 8.639 0.03; 7 11.362 0.02; 8 12.573 0.03; 9 15.691 0.05; 10 16.806 0.04; 11 20.102 0.07; 12 21.278 0.05; 13 24.787 0.09; 14 25.923 0.05; 15 29.593 0.1; 16 30.876 0.06; 17 34.634 0.12; 18 35.971 0.07; 19 39.820 0.15; 20 41.302 0.08; 21 45.474 0.15; 22 46.889 0.09; 23 51.010 0.18; 24 52.624 0.20; 25 56.846 0.22; 26 58.545 0.18; 27 63.238 0.22; 28 64.388 0.31; 29 69.126 0.31; 30 70.927 0.3]; np_1 = doerte(:,1); en_1 = doerte(:,2); den_1 = doerte(:,3); %-------------------------------------------------------------------------- F1 = figure(1); plot(Nn,En,'bo-') hold on plot(np_0,en_0,'r^-') plot(np_1,en_1,'ks-') legend('SLDA','GFMC','FN-DMC','Location','SouthEast') xlabel('N') ylabel('E(N)') % plot(Nn,R2,'bs-') % plot(np_v,en_v,'rs') % legend('SLDA E','SLDA ','GFMC E','FN-DMC','GFMC ',... % 'Location','SouthEast') hold off D_En_0 = en_0./En(1:22)-1; D_En_1 = en_1'./En-1; axes('position',[0.25,0.6,0.3,0.35]) plot(Nn(1:22),D_En_0,'r^-',np_1,D_En_1,'ks-',[0,30],zeros(2,1),'k') xlabel('N') ylabel('\delta E(N)') F2 = figure(2); Gap_s = (En(3:2:29) - (En(2:2:28)+En(4:2:30))/2)'; Gap_0 = (en_0(3:2:21) - (en_0(2:2:20)+en_0(4:2:22))/2)'; Gap_1 = (en_1(3:2:29) - (en_1(2:2:28)+en_1(4:2:30))/2); d_Gap_0 = (den_0(3:2:21)+(den_0(2:2:20)+den_0(4:2:22))/2)'; d_Gap_1 = (den_1(3:2:29)+(den_1(2:2:28)+den_1(4:2:30))/2); plot(Nn(3:2:29)',Gap_s,'bo-') hold on errorbar(Nn(3:2:21)',Gap_0,d_Gap_0,'r^-') errorbar(Nn(3:2:29)',Gap_1,d_Gap_1,'ks-') xlabel('N') ylabel('\Delta(N)') legend('SLDA','GFMC','FN-DMC','Location','NorthWest') hold off