112 lines
3.5 KiB
Matlab
112 lines
3.5 KiB
Matlab
function [psi] = runGradientDescent(this,psi,Params,Transf,VDk,V,Observ)
|
|
|
|
format long;
|
|
|
|
% Convergence Criteria:
|
|
epsilon = 1E-8;
|
|
alpha = 1E-4;
|
|
beta = 0.9;
|
|
|
|
Observ.residual = 1;
|
|
psi_old = psi; % Previous psi value (for heavy-ball method)
|
|
|
|
if this.PlotLive
|
|
Plotter.plotLive(psi,Params,Transf,Observ)
|
|
drawnow
|
|
end
|
|
|
|
% Minimization Loop
|
|
for idx = 1:this.MaxIterationsForGD
|
|
|
|
% Compute gradient
|
|
J = compute_gradient(psi, Params, Transf, VDk, V);
|
|
|
|
% Calculate chemical potential and norm
|
|
muchem = sum(real(conj(psi) .* J)) / sum(abs(psi).^2);
|
|
|
|
% Calculate residual and check convergence
|
|
residual = sum(abs(J - muchem * psi).^2) * Transf.dx * Transf.dy * Transf.dz;
|
|
if residual < epsilon
|
|
fprintf('Convergence reached at iteration %d\n', idx);
|
|
break;
|
|
end
|
|
|
|
% Update psi using heavy-ball method
|
|
psi_new = (1 + beta) * psi - alpha * J - beta * psi_old;
|
|
psi_old = psi;
|
|
|
|
% Normalize psi
|
|
Norm = sum(abs(psi_new(:)).^2) * Transf.dx * Transf.dy * Transf.dz;
|
|
psi = sqrt(Params.N) * psi_new / sqrt(Norm);
|
|
|
|
if mod(idx,500) == 0
|
|
|
|
% Collect change in energy
|
|
E = this.Calculator.calculateTotalEnergy(psi,Params,Transf,VDk,V);
|
|
E = E/Norm;
|
|
Observ.EVec = [Observ.EVec E];
|
|
|
|
% Collect chemical potentials
|
|
Observ.mucVec = [Observ.mucVec muchem];
|
|
|
|
% Collect residuals
|
|
Observ.residual = [Observ.residual res];
|
|
Observ.res_idx = Observ.res_idx + 1;
|
|
|
|
if this.PlotLive
|
|
Plotter.plotLive(psi,Params,Transf,Observ)
|
|
drawnow
|
|
end
|
|
|
|
save(sprintf(strcat(this.SaveDirectory, '/Run_%03i/psi_gs.mat'),Params.njob),'psi','muchem','Observ','Transf','Params','VDk','V');
|
|
end
|
|
end
|
|
|
|
% Check if max iterations were hit without convergence
|
|
last_iteration_number = idx;
|
|
if last_iteration_number == this.MaxIterationsForGD
|
|
fprintf('Max iterations reached without convergence. Final chemical potential: %.6f, Residual: %.6f\n', muchem, residual);
|
|
else
|
|
fprintf('Converged in %d iterations. Final chemical potential: %.6f\n', last_iteration_number, muchem);
|
|
end
|
|
|
|
% Change in Energy
|
|
E = this.Calculator.calculateTotalEnergy(psi,Params,Transf,VDk,V);
|
|
E = E/Norm;
|
|
Observ.EVec = [Observ.EVec E];
|
|
|
|
disp('Saving data...');
|
|
save(sprintf(strcat(this.SaveDirectory, '/Run_%03i/psi_gs.mat'),Params.njob),'psi','muchem','Observ','Transf','Params','VDk','V');
|
|
disp('Save complete!');
|
|
|
|
end
|
|
|
|
%% Helper functions
|
|
% Numerical Gradient Calculation using the finite differences method
|
|
function J = compute_gradient(psi, Params, Transf, VDk, V)
|
|
|
|
% Operators
|
|
|
|
% Kinetic energy
|
|
KEop = 0.5 * (Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
|
|
HKin = @(w) real(ifft(KEop.*fft(w)));
|
|
|
|
% Trap Potential
|
|
HV = @(w) V.*w;
|
|
|
|
% Contact interactions
|
|
Hint = @(w) (Params.gs*abs(psi).^2).*w;
|
|
|
|
% DDIs
|
|
frho = fftn(abs(psi).^2);
|
|
Phi = real(ifftn(frho.*VDk));
|
|
Hddi = @(w) (Params.gdd*Phi).*w;
|
|
|
|
% Quantum fluctuations
|
|
Hqf = @(w) (Params.gammaQF*abs(psi).^3).*w;
|
|
|
|
H = @(w) HKin(w) + HV(w) + Hint(w) + Hddi(w) + Hqf(w);
|
|
|
|
J = H(psi);
|
|
|
|
end |