Restored individual files for functions of Calculator class.
This commit is contained in:
parent
1e94302160
commit
3a4419fe67
@ -0,0 +1,29 @@
|
|||||||
|
function muchem = calculateChemicalPotential(~,psi,Params,Transf,VDk,V)
|
||||||
|
|
||||||
|
%Parameters
|
||||||
|
normfac = Params.Lx*Params.Ly*Params.Lz/numel(psi);
|
||||||
|
KEop= 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
|
||||||
|
|
||||||
|
% DDIs
|
||||||
|
frho=fftn(abs(psi).^2);
|
||||||
|
Phi=real(ifftn(frho.*VDk));
|
||||||
|
|
||||||
|
Eddi = (Params.gdd*Phi.*abs(psi).^2);
|
||||||
|
|
||||||
|
%Kinetic energy
|
||||||
|
Ekin = KEop.*abs(fftn(psi)*normfac).^2;
|
||||||
|
Ekin = trapz(Ekin(:))*Transf.dkx*Transf.dky*Transf.dkz/(2*pi)^3;
|
||||||
|
|
||||||
|
%Potential energy
|
||||||
|
Epot = V.*abs(psi).^2;
|
||||||
|
|
||||||
|
%Contact interactions
|
||||||
|
Eint = Params.gs*abs(psi).^4;
|
||||||
|
|
||||||
|
%Quantum fluctuations
|
||||||
|
Eqf = Params.gammaQF*abs(psi).^5;
|
||||||
|
|
||||||
|
%Total energy
|
||||||
|
muchem = Ekin + trapz(Epot(:) + Eint(:) + Eddi(:) + Eqf(:))*Transf.dx*Transf.dy*Transf.dz; %
|
||||||
|
muchem = muchem / Params.N;
|
||||||
|
end
|
||||||
@ -0,0 +1,35 @@
|
|||||||
|
function E = calculateEnergyComponents(~,psi,Params,Transf,VDk,V)
|
||||||
|
|
||||||
|
%Parameters
|
||||||
|
|
||||||
|
KEop= 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
|
||||||
|
normfac = Params.Lx*Params.Ly*Params.Lz/numel(psi);
|
||||||
|
|
||||||
|
% DDIs
|
||||||
|
frho = fftn(abs(psi).^2);
|
||||||
|
Phi = real(ifftn(frho.*VDk));
|
||||||
|
|
||||||
|
Eddi = 0.5*Params.gdd*Phi.*abs(psi).^2;
|
||||||
|
E.Eddi = trapz(Eddi(:))*Transf.dx*Transf.dy*Transf.dz;
|
||||||
|
|
||||||
|
% EddiTot = trapz(Eddi(:))*Transf.dx*Transf.dy*Transf.dz;
|
||||||
|
|
||||||
|
%Kinetic energy
|
||||||
|
% psik = ifftshift(fftn(fftshift(psi)))*normfac;
|
||||||
|
|
||||||
|
Ekin = KEop.*abs(fftn(psi)*normfac).^2;
|
||||||
|
E.Ekin = trapz(Ekin(:))*Transf.dkx*Transf.dky*Transf.dkz/(2*pi)^3;
|
||||||
|
|
||||||
|
% Potential energy
|
||||||
|
Epot = V.*abs(psi).^2;
|
||||||
|
E.Epot = trapz(Epot(:))*Transf.dx*Transf.dy*Transf.dz;
|
||||||
|
|
||||||
|
%Contact interactions
|
||||||
|
Eint = 0.5*Params.gs*abs(psi).^4;
|
||||||
|
E.Eint = trapz(Eint(:))*Transf.dx*Transf.dy*Transf.dz;
|
||||||
|
|
||||||
|
%Quantum fluctuations
|
||||||
|
Eqf = 0.4*Params.gammaQF*abs(psi).^5;
|
||||||
|
E.Eqf = trapz(Eqf(:))*Transf.dx*Transf.dy*Transf.dz;
|
||||||
|
|
||||||
|
end
|
||||||
@ -0,0 +1,25 @@
|
|||||||
|
function res = calculateNormalizedResiduals(~,psi,Params,Transf,VDk,V,muchem)
|
||||||
|
|
||||||
|
KEop= 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
|
||||||
|
|
||||||
|
% DDIs
|
||||||
|
frho=fftn(abs(psi).^2);
|
||||||
|
Phi=real(ifftn(frho.*VDk));
|
||||||
|
|
||||||
|
Eddi = Params.gdd*Phi.*psi;
|
||||||
|
|
||||||
|
%Kinetic energy
|
||||||
|
Ekin = ifftn(KEop.*fftn(psi));
|
||||||
|
|
||||||
|
%Potential energy
|
||||||
|
Epot = V.*psi;
|
||||||
|
|
||||||
|
%Contact interactions
|
||||||
|
Eint = Params.gs*abs(psi).^2.*psi;
|
||||||
|
|
||||||
|
%Quantum fluctuations
|
||||||
|
Eqf = Params.gammaQF*abs(psi).^3.*psi;
|
||||||
|
|
||||||
|
%Total energy
|
||||||
|
res = trapz(abs(Ekin(:) + Epot(:) + Eint(:) + Eddi(:) + Eqf(:) - muchem*psi(:))*Transf.dx*Transf.dy*Transf.dz)/trapz(abs(muchem*psi(:))*Transf.dx*Transf.dy*Transf.dz);
|
||||||
|
end
|
||||||
@ -0,0 +1,39 @@
|
|||||||
|
function VDkSemi = calculateNumericalHankelTransform(~,kr,kz,Rmax,Zmax,Nr)
|
||||||
|
|
||||||
|
% accuracy inputs for numerical integration
|
||||||
|
if(nargin==5)
|
||||||
|
Nr = 5e4;
|
||||||
|
end
|
||||||
|
|
||||||
|
Nz = 64;
|
||||||
|
farRmultiple = 2000;
|
||||||
|
|
||||||
|
% midpoint grids for the integration over 0<z<Zmax, Rmax<r<farRmultiple*Rmax (i.e. starts at Rmax)
|
||||||
|
dr=(farRmultiple-1)*Rmax/Nr;
|
||||||
|
r = ((1:Nr)'-0.5)*dr+Rmax;
|
||||||
|
dz=Zmax/Nz;
|
||||||
|
z = ((1:Nz)-0.5)*dz;
|
||||||
|
[R, Z] = ndgrid(r,z);
|
||||||
|
Rsq = R.^2 + Z.^2;
|
||||||
|
|
||||||
|
% real space interaction to be transformed
|
||||||
|
igrandbase = (1 - 3*Z.^2./Rsq)./Rsq.^(3/2);
|
||||||
|
|
||||||
|
% do the Hankel/Fourier-Bessel transform numerically
|
||||||
|
% prestore to ensure each besselj is only calculated once
|
||||||
|
% cell is faster than (:,:,krn) slicing
|
||||||
|
Nkr = numel(kr);
|
||||||
|
besselr = cell(Nkr,1);
|
||||||
|
for krn = 1:Nkr
|
||||||
|
besselr{krn} = repmat(r.*besselj(0,kr(krn)*r),1,Nz);
|
||||||
|
end
|
||||||
|
|
||||||
|
VDkSemi = zeros([Nkr,numel(kz)]);
|
||||||
|
for kzn = 1:numel(kz)
|
||||||
|
igrandbasez = repmat(cos(kz(kzn)*z),Nr,1) .* igrandbase;
|
||||||
|
for krn = 1:Nkr
|
||||||
|
igrand = igrandbasez.*besselr{krn};
|
||||||
|
VDkSemi(krn,kzn) = VDkSemi(krn,kzn) - sum(igrand(:))*dz*dr;
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
@ -0,0 +1,58 @@
|
|||||||
|
function [m_Order] = calculateOrderParameter(~,psi,Transf,Params,VDk,V,T,muchem)
|
||||||
|
|
||||||
|
NumRealiz = 100;
|
||||||
|
|
||||||
|
Mx = numel(Transf.x);
|
||||||
|
My = numel(Transf.y);
|
||||||
|
Mz = numel(Transf.z);
|
||||||
|
|
||||||
|
r = normrnd(0,1,size(psi));
|
||||||
|
theta = rand(size(psi));
|
||||||
|
noise = r.*exp(2*pi*1i*theta);
|
||||||
|
|
||||||
|
KEop= 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
|
||||||
|
Gamma = 1-1i*Params.gamma_S;
|
||||||
|
dt = Params.dt;
|
||||||
|
|
||||||
|
avgpsi = 0;
|
||||||
|
avgpsi2 = 0;
|
||||||
|
|
||||||
|
for jj = 1:NumRealiz
|
||||||
|
%generate initial state
|
||||||
|
xi = sqrt(2*Params.gamma_S*Params.kbol*T*10^(-9)*dt/(Params.hbar*Params.w0*Transf.dx*Transf.dy*Transf.dz));
|
||||||
|
swapx = randi(length(Transf.x),1,length(Transf.x));
|
||||||
|
swapy = randi(length(Transf.y),1,length(Transf.y));
|
||||||
|
swapz = randi(length(Transf.z),1,length(Transf.z));
|
||||||
|
psi_j = psi + xi * noise(swapx,swapy,swapz);
|
||||||
|
|
||||||
|
% --- % propagate forward in time 1 time step:
|
||||||
|
%kin
|
||||||
|
psi_j = fftn(psi_j);
|
||||||
|
psi_j = psi_j.*exp(-0.5*1i*Gamma*dt*KEop);
|
||||||
|
psi_j = ifftn(psi_j);
|
||||||
|
|
||||||
|
%DDI
|
||||||
|
frho = fftn(abs(psi_j).^2);
|
||||||
|
Phi = real(ifftn(frho.*VDk));
|
||||||
|
|
||||||
|
%Real-space
|
||||||
|
psi_j = psi_j.*exp(-1i*Gamma*dt*(V + Params.gs*abs(psi_j).^2 + Params.gammaQF*abs(psi_j).^3 + Params.gdd*Phi - muchem));
|
||||||
|
|
||||||
|
%kin
|
||||||
|
psi_j = fftn(psi_j);
|
||||||
|
psi_j = psi_j.*exp(-0.5*1i*Gamma*dt*KEop);
|
||||||
|
psi_j = ifftn(psi_j);
|
||||||
|
|
||||||
|
%Projection
|
||||||
|
kcut = sqrt(2*Params.e_cut);
|
||||||
|
K = (Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2)<kcut.^2;
|
||||||
|
psi_j = ifftn(K.*fftn(psi_j));
|
||||||
|
|
||||||
|
% --- %
|
||||||
|
avgpsi = avgpsi + abs(sum(psi_j(:)))/NumRealiz;
|
||||||
|
avgpsi2 = avgpsi2 + sum(abs(psi_j(:)).^2)/NumRealiz;
|
||||||
|
|
||||||
|
end
|
||||||
|
|
||||||
|
m_Order = 1/sqrt(Mx*My*Mz)*avgpsi/sqrt(avgpsi2);
|
||||||
|
end
|
||||||
@ -0,0 +1,19 @@
|
|||||||
|
function [PhaseC] = calculatePhaseCoherence(~,psi,Transf,Params)
|
||||||
|
|
||||||
|
norm = sum(sum(sum(abs(psi).^2,1),2),3)*Transf.dx*Transf.dy*Transf.dz;
|
||||||
|
psi = psi/sqrt(norm);
|
||||||
|
|
||||||
|
NumGlobalShifts = 800;
|
||||||
|
betas = []; phishift = [];
|
||||||
|
for jj = 1:NumGlobalShifts
|
||||||
|
phishift(jj) = -pi + 2*pi*(jj-1)/NumGlobalShifts;
|
||||||
|
betas(jj) = sum(sum(sum(abs(angle(psi*exp(-1i*phishift(jj)))).*abs(psi).^2)));
|
||||||
|
end
|
||||||
|
[minbeta,minidx] = min(betas);
|
||||||
|
|
||||||
|
psi = psi*exp(-1i*phishift(minidx));
|
||||||
|
phi = angle(psi);
|
||||||
|
|
||||||
|
avgphi = sum(sum(sum(phi.*abs(psi).^2,1),2),3)*Transf.dx*Transf.dy*Transf.dz;
|
||||||
|
PhaseC = sum(sum(sum(abs(angle(psi)-avgphi).*abs(psi).^2)))*Transf.dx*Transf.dy*Transf.dz;
|
||||||
|
end
|
||||||
@ -0,0 +1,32 @@
|
|||||||
|
function E = calculateTotalEnergy(~,psi,Params,Transf,VDk,V)
|
||||||
|
|
||||||
|
%Parameters
|
||||||
|
|
||||||
|
KEop= 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
|
||||||
|
normfac = Params.Lx*Params.Ly*Params.Lz/numel(psi);
|
||||||
|
|
||||||
|
% DDIs
|
||||||
|
frho = fftn(abs(psi).^2);
|
||||||
|
Phi = real(ifftn(frho.*VDk));
|
||||||
|
|
||||||
|
Eddi = 0.5*Params.gdd*Phi.*abs(psi).^2;
|
||||||
|
|
||||||
|
% EddiTot = trapz(Eddi(:))*Transf.dx*Transf.dy*Transf.dz;
|
||||||
|
|
||||||
|
%Kinetic energy
|
||||||
|
% psik = ifftshift(fftn(fftshift(psi)))*normfac;
|
||||||
|
|
||||||
|
Ekin = KEop.*abs(fftn(psi)*normfac).^2;
|
||||||
|
Ekin = trapz(Ekin(:))*Transf.dkx*Transf.dky*Transf.dkz/(2*pi)^3;
|
||||||
|
|
||||||
|
% Potential energy
|
||||||
|
Epot = V.*abs(psi).^2;
|
||||||
|
|
||||||
|
%Contact interactions
|
||||||
|
Eint = 0.5*Params.gs*abs(psi).^4;
|
||||||
|
|
||||||
|
%Quantum fluctuations
|
||||||
|
Eqf = 0.4*Params.gammaQF*abs(psi).^5;
|
||||||
|
|
||||||
|
E = Ekin + trapz(Epot(:) + Eint(:) + Eddi(:) + Eqf(:))*Transf.dx*Transf.dy*Transf.dz;
|
||||||
|
end
|
||||||
@ -0,0 +1,64 @@
|
|||||||
|
function VDk = calculateVDCutoff(this,Params,Transf,TransfRad)
|
||||||
|
% makes the dipolar interaction matrix, size numel(Params.kr) * numel(Params.kz)
|
||||||
|
% Rmax and Zmax are the interaction cutoffs
|
||||||
|
% VDk needs to be multiplied by Cdd
|
||||||
|
% approach is that of Lu, PRA 82, 023622 (2010)
|
||||||
|
|
||||||
|
% == Calulating the DDI potential in Fourier space with appropriate cutoff == %
|
||||||
|
% Cylindrical (semianalytic)
|
||||||
|
% Cylindrical infinite Z, polarized along x (analytic)
|
||||||
|
% Spherical
|
||||||
|
|
||||||
|
switch this.CutoffType
|
||||||
|
case 'Cylindrical' %Cylindrical (semianalytic)
|
||||||
|
Zcutoff = Params.Lz/2;
|
||||||
|
alph = acos((Transf.KX*sin(Params.theta)*cos(Params.phi)+Transf.KY*sin(Params.theta)*sin(Params.phi)+Transf.KZ*cos(Params.theta))./sqrt(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2));
|
||||||
|
alph(1) = pi/2;
|
||||||
|
|
||||||
|
% Analytic part of cutoff for slice 0<z<Zmax, 0<r<Inf Ronen, PRL 98, 030406 (2007)
|
||||||
|
cossq = cos(alph).^2;
|
||||||
|
VDk = cossq-1/3;
|
||||||
|
sinsq = 1 - cossq;
|
||||||
|
VDk = VDk + exp(-Zcutoff*sqrt(Transf.KX.^2+Transf.KY.^2)).*( sinsq .* cos(Zcutoff * Transf.KZ) - sqrt(sinsq.*cossq).*sin(Zcutoff * Transf.KZ) );
|
||||||
|
|
||||||
|
% Nonanalytic part
|
||||||
|
% For a cylindrical cutoff, we need to construct a kr grid based on the 3D parameters using Bessel quadrature
|
||||||
|
VDkNon = this.calculateNumericalHankelTransform(TransfRad.kr, TransfRad.kz, TransfRad.Rmax, Zcutoff);
|
||||||
|
|
||||||
|
% Interpolating the nonanalytic part onto 3D grid
|
||||||
|
fullkr = [-flip(TransfRad.kr)',TransfRad.kr'];
|
||||||
|
[KR,KZ] = ndgrid(fullkr,TransfRad.kz);
|
||||||
|
[KX3D,KY3D,KZ3D] = ndgrid(ifftshift(Transf.kx),ifftshift(Transf.ky),ifftshift(Transf.kz));
|
||||||
|
KR3D = sqrt(KX3D.^2 + KY3D.^2);
|
||||||
|
fullVDK = [flip(VDkNon',2),VDkNon']';
|
||||||
|
VDkNon = interpn(KR,KZ,fullVDK,KR3D,KZ3D,'spline',0); %Last argument is -1/3 for full VDk. 0 for nonanalytic piece?
|
||||||
|
VDkNon = fftshift(VDkNon);
|
||||||
|
|
||||||
|
VDk = VDk + VDkNon;
|
||||||
|
|
||||||
|
case 'CylindricalInfiniteZ' %Cylindrical infinite Z, polarized along x -- PRA 107, 033301 (2023)
|
||||||
|
alph = acos((Transf.KX*sin(Params.theta)*cos(Params.phi)+Transf.KY*sin(Params.theta)*sin(Params.phi)+Transf.KZ*cos(Params.theta))./sqrt(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2));
|
||||||
|
alph(1) = pi/2;
|
||||||
|
rhoc = max([abs(Transf.x),abs(Transf.y)]);
|
||||||
|
KR = sqrt(Transf.KX.^2+Transf.KY.^2);
|
||||||
|
func = @(n,u,v) v.^2./(u.^2+v.^2).*(v.*besselj(n,u).*besselk(n+1,v) - u.*besselj(n+1,u).*besselk(n,v));
|
||||||
|
VDk = -0.5*func(0,KR*rhoc,abs(Transf.KZ)*rhoc) + (Transf.KX.^2./KR.^2 - 0.5).*func(2,KR*rhoc,abs(Transf.KZ)*rhoc);
|
||||||
|
VDk = (1/3)*(3*(cos(alph).^2)-1) - VDk;
|
||||||
|
|
||||||
|
VDk(KR==0) = -1/3 + 1/2*abs(Transf.KZ(KR==0))*rhoc.*besselk(1,abs(Transf.KZ(KR==0))*rhoc);
|
||||||
|
VDk(Transf.KZ==0) = 1/6 + (Transf.KX(Transf.KZ==0).^2-Transf.KY(Transf.KZ==0).^2)./(KR(Transf.KZ==0).^2).*(1/2 - besselj(1,KR(Transf.KZ==0)*rhoc)./(KR(Transf.KZ==0)*rhoc));
|
||||||
|
VDk(1,1,1) = 1/6;
|
||||||
|
|
||||||
|
case 'Spherical' %Spherical
|
||||||
|
Rcut = min(Params.Lx/2,Params.Ly/2,Params.Lz/2);
|
||||||
|
alph = acos((Transf.KX*sin(Params.theta)*cos(Params.phi)+Transf.KY*sin(Params.theta)*sin(Params.phi)+Transf.KZ*cos(Params.theta))./sqrt(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2));
|
||||||
|
alph(1) = pi/2;
|
||||||
|
|
||||||
|
K = sqrt(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
|
||||||
|
VDk = (cos(alph).^2-1/3).*(1 + 3*cos(Rcut*K)./(Rcut^2.*K.^2) - 3*sin(Rcut*K)./(Rcut^3.*K.^3));
|
||||||
|
|
||||||
|
otherwise
|
||||||
|
disp('Choose a valid DDI cutoff type!')
|
||||||
|
return
|
||||||
|
end
|
||||||
|
end
|
||||||
Loading…
Reference in New Issue
Block a user