218 lines
8.7 KiB
Matlab
218 lines
8.7 KiB
Matlab
function [psi] = propagateWavefunction(this,psi,Params,Transf,VDk,V,t_idx,Observ)
|
|
|
|
switch this.SimulationMode
|
|
case 'ImaginaryTimeEvolution'
|
|
|
|
dt = -1j*abs(Params.dt); % Imaginary Time
|
|
|
|
KEop = 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
|
|
Observ.residual = 1;
|
|
Observ.res = 1;
|
|
|
|
muchem = this.Calculator.calculateChemicalPotential(psi,Params,Transf,VDk,V);
|
|
|
|
if this.PlotLive
|
|
Plotter.plotLive(psi,Params,Transf,Observ,this.SimulationMode)
|
|
drawnow
|
|
end
|
|
|
|
AdaptIdx = 0;
|
|
|
|
while t_idx < Params.sim_time_cut_off
|
|
|
|
% kin
|
|
psi = fftn(psi);
|
|
psi = psi.*exp(-0.5*1i*dt*KEop);
|
|
psi = ifftn(psi);
|
|
|
|
% DDI
|
|
frho = fftn(abs(psi).^2);
|
|
Phi = real(ifftn(frho.*VDk));
|
|
|
|
% Real-space
|
|
psi = psi.*exp(-1i*dt*(V + Params.gs*abs(psi).^2 + Params.gammaQF*abs(psi).^3 + Params.gdd*Phi - muchem));
|
|
|
|
% kin
|
|
psi = fftn(psi);
|
|
psi = psi.*exp(-0.5*1i*dt*KEop);
|
|
psi = ifftn(psi);
|
|
|
|
% Renorm
|
|
Norm = sum(abs(psi(:)).^2)*Transf.dx*Transf.dy*Transf.dz;
|
|
psi = sqrt(Params.N)*psi/sqrt(Norm);
|
|
|
|
muchem = this.Calculator.calculateChemicalPotential(psi,Params,Transf,VDk,V);
|
|
|
|
if mod(t_idx,1000) == 0
|
|
|
|
% Change in Energy
|
|
E = this.Calculator.calculateTotalEnergy(psi,Params,Transf,VDk,V);
|
|
E = E/Norm;
|
|
Observ.EVec = [Observ.EVec E];
|
|
|
|
% Chemical potential
|
|
Observ.mucVec = [Observ.mucVec muchem];
|
|
|
|
% Normalized residuals
|
|
res = this.Calculator.calculateNormalizedResiduals(psi,Params,Transf,VDk,V,muchem);
|
|
Observ.residual = [Observ.residual res];
|
|
Observ.res_idx = Observ.res_idx + 1;
|
|
|
|
if this.PlotLive
|
|
Plotter.plotLive(psi,Params,Transf,Observ,this.SimulationMode)
|
|
drawnow
|
|
end
|
|
|
|
if this.DoSave
|
|
save(sprintf(strcat(this.SaveDirectory, '/Run_%03i/psi_gs.mat'),Params.njob),'psi','muchem','Observ','t_idx','Transf','Params','VDk','V');
|
|
end
|
|
|
|
%Adaptive time step -- Careful, this can quickly get out of control
|
|
relres = abs(Observ.residual(Observ.res_idx)-Observ.residual(Observ.res_idx-1))/Observ.residual(Observ.res_idx);
|
|
if relres < 1e-5
|
|
if AdaptIdx > 4 && abs(dt) > Params.mindt
|
|
dt = dt / 2;
|
|
fprintf('Time step changed to '); disp(dt);
|
|
AdaptIdx = 0;
|
|
elseif AdaptIdx > 4 && abs(dt) < Params.mindt
|
|
break
|
|
else
|
|
AdaptIdx = AdaptIdx + 1;
|
|
end
|
|
else
|
|
AdaptIdx = 0;
|
|
end
|
|
end
|
|
if any(isnan(psi(:)))
|
|
disp('NaNs encountered!')
|
|
break
|
|
end
|
|
t_idx=t_idx+1;
|
|
end
|
|
|
|
% Change in Energy
|
|
E = this.Calculator.calculateTotalEnergy(psi,Params,Transf,VDk,V);
|
|
E = E/Norm;
|
|
Observ.EVec = [Observ.EVec E];
|
|
|
|
Observ.res_idx = Observ.res_idx + 1;
|
|
|
|
disp('Run completed!');
|
|
if this.DoSave
|
|
disp('Saving data...');
|
|
save(sprintf(strcat(this.SaveDirectory, '/Run_%03i/psi_gs.mat'),Params.njob),'psi','muchem','Observ','t_idx','Transf','Params','VDk','V');
|
|
disp('Save complete!');
|
|
end
|
|
|
|
case 'RealTimeEvolution'
|
|
|
|
dt = abs(Params.dt);
|
|
|
|
KEop = 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
|
|
|
|
muchem = this.Calculator.calculateChemicalPotential(psi,Params,Transf,VDk,V);
|
|
|
|
if this.PlotLive
|
|
Plotter.plotLive(psi,Params,Transf,Observ,this.SimulationMode)
|
|
drawnow
|
|
end
|
|
|
|
while t_idx < this.QuenchSettings.tSteps
|
|
|
|
%Time
|
|
tVal = this.QuenchSettings.tVec(t_idx);
|
|
|
|
if this.QuenchScatteringLength
|
|
Params.gs = this.QuenchSettings.gs_vec(t_idx);
|
|
Params.gammaQF = this.QuenchSettings.gammaQF_vec(t_idx);
|
|
end
|
|
|
|
if this.RotateDipoles
|
|
Params.theta = this.QuenchSettings.theta_vec(t_idx);
|
|
Params.phi = this.QuenchSettings.phi_vec(t_idx);
|
|
VDk = this.Calculator.calculateVDk(Params, Transf, TransfRad, this.IncludeDDICutOff);
|
|
end
|
|
|
|
% Parameters at time t
|
|
|
|
%kin
|
|
psi = fftn(psi);
|
|
psi = psi.*exp(-0.5*1i*dt*KEop);
|
|
psi = ifftn(psi);
|
|
|
|
%DDI
|
|
frho = fftn(abs(psi).^2);
|
|
Phi = real(ifftn(frho.*VDk));
|
|
|
|
%Real-space
|
|
psi = psi.*exp(-1i*dt*(V + Params.gs*abs(psi).^2 + Params.gammaQF*abs(psi).^3 + Params.gdd*Phi - muchem));
|
|
|
|
%kin
|
|
psi = fftn(psi);
|
|
psi = psi.*exp(-0.5*1i*dt*KEop);
|
|
psi = ifftn(psi);
|
|
|
|
muchem = this.Calculator.calculateChemicalPotential(psi,Params,Transf,VDk,V);
|
|
|
|
if mod(t_idx,this.QuenchSettings.saveinterval)==0
|
|
%Change in Normalization
|
|
Norm = sum(abs(psi(:)).^2)*Transf.dx*Transf.dy*Transf.dz; % normalisation
|
|
Observ.NormVec = [Observ.NormVec Norm];
|
|
|
|
%Change in Energy
|
|
E = this.Calculator.calculateTotalEnergy(psi,Params,Transf,VDk,V);
|
|
E = E/Norm;
|
|
Observ.EVec = [Observ.EVec E];
|
|
|
|
% Phase coherence
|
|
[PhaseC] = this.Calculator.calculatePhaseCoherence(psi,Transf);
|
|
Observ.PCVec = [Observ.PCVec PhaseC];
|
|
|
|
Observ.tVecPlot = [Observ.tVecPlot tVal];
|
|
Observ.res_idx = Observ.res_idx + 1;
|
|
|
|
if this.PlotLive
|
|
Plotter.plotLive(psi,Params,Transf,Observ,this.SimulationMode)
|
|
drawnow
|
|
end
|
|
|
|
if this.DoSave
|
|
save(sprintf(strcat(this.SaveDirectory, '/Run_%03i/psi_%i.mat'), Params.njob, Observ.res_idx),'psi','muchem','Observ','t_idx','Transf','Params');
|
|
end
|
|
end
|
|
if any(isnan(psi(:)))
|
|
disp('NaNs encountered!')
|
|
break
|
|
end
|
|
t_idx = t_idx+1;
|
|
end
|
|
|
|
%Change in Normalization
|
|
Norm = sum(abs(psi(:)).^2)*Transf.dx*Transf.dy*Transf.dz; % normalisation
|
|
Observ.NormVec = [Observ.NormVec Norm];
|
|
|
|
%Change in Energy
|
|
E = this.Calculator.calculateTotalEnergy(psi,Params,Transf,VDk,V);
|
|
E = E/Norm;
|
|
Observ.EVec = [Observ.EVec E];
|
|
|
|
% Phase coherence
|
|
[PhaseC] = this.Calculator.calculatePhaseCoherence(psi,Transf);
|
|
Observ.PCVec = [Observ.PCVec PhaseC];
|
|
|
|
Observ.tVecPlot = [Observ.tVecPlot tVal];
|
|
Observ.res_idx = Observ.res_idx + 1;
|
|
|
|
disp('Run completed!');
|
|
if this.DoSave
|
|
disp('Saving data...');
|
|
save(sprintf(strcat(this.SaveDirectory, '/Run_%03i/psi_%i.mat'),Params.njob,Observ.res_idx),'psi','muchem','Observ','t_idx','Transf','Params');
|
|
disp('Save complete!');
|
|
end
|
|
|
|
|
|
otherwise
|
|
disp('Choose a valid method!')
|
|
return
|
|
end
|
|
end |