293 lines
12 KiB
Matlab
293 lines
12 KiB
Matlab
function [psi] = runGradientDescent(this,psi,Params,Transf,VDk,V,Observ)
|
|
|
|
format long;
|
|
|
|
AdditionalParams.GradientDescentMethod = this.GradientDescentMethod;
|
|
SimulationMode = this.SimulationMode;
|
|
|
|
switch this.GradientDescentMethod
|
|
case 'HeavyBall'
|
|
% Convergence Criteria:
|
|
alpha = 1E-5;
|
|
beta = 0.9;
|
|
epsilon = 1E-6;
|
|
|
|
Observ.residual = 1;
|
|
Observ.res = 1;
|
|
psi_old = psi; % Previous psi value (for heavy-ball method)
|
|
|
|
% Live Plotter
|
|
if this.PlotLive
|
|
Plotter.plotLive(psi,Params,Transf,Observ,SimulationMode,AdditionalParams)
|
|
drawnow
|
|
end
|
|
|
|
% Minimization Loop
|
|
for idx = 1:this.MaxIterationsForGD
|
|
|
|
% Compute gradient
|
|
J = compute_gradient(psi, Params, Transf, VDk, V);
|
|
|
|
% Calculate chemical potential
|
|
muchem = real(inner_product(J, psi, Transf)) / inner_product(psi, psi, Transf);
|
|
|
|
% Calculate residual and check convergence
|
|
diff = J - (muchem * psi);
|
|
residual = inner_product(diff, diff, Transf);
|
|
|
|
if residual < epsilon
|
|
fprintf('Convergence reached at iteration %d\n', idx);
|
|
break;
|
|
else
|
|
|
|
% Update psi using heavy-ball method
|
|
psi_new = (1 + beta) * psi - alpha * J - beta * psi_old;
|
|
psi_old = psi;
|
|
|
|
% Normalize psi
|
|
Norm = inner_product(psi_new, psi_new, Transf);
|
|
psi = psi_new / sqrt(Norm);
|
|
|
|
% Write output at specified intervals
|
|
if mod(idx,100) == 0
|
|
|
|
% Collect change in energy
|
|
E = this.Calculator.calculateTotalEnergy(psi,Params,Transf,VDk,V) / Norm;
|
|
Observ.EVec = [Observ.EVec E];
|
|
|
|
% Collect chemical potentials
|
|
Observ.mucVec = [Observ.mucVec muchem];
|
|
|
|
% Collect 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,SimulationMode,AdditionalParams)
|
|
drawnow
|
|
end
|
|
|
|
save(sprintf(strcat(this.SaveDirectory, '/Run_%03i/psi_gs.mat'),Params.njob),'psi','muchem','Observ','Transf','Params','VDk','V');
|
|
end
|
|
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) / Norm;
|
|
Observ.EVec = [Observ.EVec E];
|
|
|
|
save(sprintf(strcat(this.SaveDirectory, '/Run_%03i/psi_gs.mat'),Params.njob),'psi','muchem','Observ','Transf','Params','VDk','V');
|
|
disp('Completed and saved!');
|
|
|
|
case 'NonLinearCGD'
|
|
|
|
% Convergence Criteria:
|
|
epsilon = 1E-13;
|
|
|
|
% Iteration Counter:
|
|
i = 1;
|
|
Observ.residual = 1;
|
|
Observ.res = 1;
|
|
Observ.theta = NaN;
|
|
|
|
% Initialize the PrematureExitFlag to false
|
|
PrematureExitFlag = false;
|
|
|
|
% Live plotter
|
|
if this.PlotLive
|
|
Plotter.plotLive(psi,Params,Transf,Observ,SimulationMode,AdditionalParams)
|
|
drawnow
|
|
end
|
|
|
|
% Minimization Loop
|
|
while true
|
|
% Compute gradient
|
|
J = compute_gradient(psi, Params, Transf, VDk, V);
|
|
% Calculate chemical potential
|
|
muchem = real(inner_product(J, psi, Transf)) / inner_product(psi, psi, Transf);
|
|
% Calculate residual
|
|
residual = J - (muchem * psi);
|
|
|
|
% Energy convergence check every 100 steps
|
|
if mod(i,100) == 0 && length(Observ.EVec) > 1
|
|
energydifference = abs(Observ.EVec(end) - Observ.EVec(end-1));
|
|
if energydifference <= epsilon
|
|
disp('Tolerance reached: Energy difference is below the specified epsilon.');
|
|
PrematureExitFlag = true;
|
|
break;
|
|
end
|
|
end
|
|
|
|
if i >= this.MaxIterationsForGD
|
|
disp('Maximum number of iterations for CGD reached.');
|
|
PrematureExitFlag = true;
|
|
break;
|
|
end
|
|
|
|
% Initialize search direction if first iteration
|
|
if i == 1
|
|
d = -residual;
|
|
else % Compute beta via Polak-Ribiere and create new direction
|
|
residual_new = residual;
|
|
beta = compute_beta(residual_new, residual_old, Transf);
|
|
d = -residual_new + beta * p_old;
|
|
end
|
|
|
|
residual_old = residual;
|
|
proj = inner_product(d, psi, Transf);
|
|
p = d - (proj * psi);
|
|
p_old = p;
|
|
|
|
% Compute optimal theta to generate psi in the direction of minimum in the energy landscape
|
|
theta = compute_optimal_theta(p, muchem, psi, Params, Transf, VDk, V);
|
|
|
|
% Update solution
|
|
gamma = 1 / sqrt(inner_product(p, p, Transf));
|
|
psi = (cos(theta).*psi) + (sin(theta).*(p*gamma));
|
|
|
|
% Normalize psi
|
|
Norm = abs(inner_product(psi, psi, Transf));
|
|
psi = psi / sqrt(Norm);
|
|
i = i + 1;
|
|
|
|
% Calculate chemical potential with new psi
|
|
J = compute_gradient(psi, Params, Transf, VDk, V);
|
|
muchem = real(inner_product(J, psi, Transf)) / Norm;
|
|
|
|
if mod(i,100) == 0
|
|
|
|
% Collect Energy value
|
|
E = this.Calculator.calculateTotalEnergy(psi,Params,Transf,VDk,V) / Norm;
|
|
Observ.EVec = [Observ.EVec E];
|
|
|
|
% Collect Chemical potential value
|
|
Observ.mucVec = [Observ.mucVec muchem];
|
|
|
|
% Collect Normalized residuals
|
|
res = this.Calculator.calculateNormalizedResiduals(psi,Params,Transf,VDk,V,muchem);
|
|
Observ.residual = [Observ.residual res];
|
|
Observ.res_idx = Observ.res_idx + 1;
|
|
|
|
% Collect calculated thetas
|
|
Observ.theta = [Observ.theta theta];
|
|
|
|
% Live plotter
|
|
if this.PlotLive
|
|
Plotter.plotLive(psi,Params,Transf,Observ,SimulationMode,AdditionalParams)
|
|
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 loop ended prematurely
|
|
if PrematureExitFlag
|
|
disp('Optimizer ended prematurely without convergence to a minimum.');
|
|
else
|
|
fprintf('Minimum found! Number of Iterations for Convergence: %d\n\n', i);
|
|
end
|
|
|
|
% Change in Energy
|
|
E = this.Calculator.calculateTotalEnergy(psi,Params,Transf,VDk,V) / Norm;
|
|
Observ.EVec = [Observ.EVec E];
|
|
|
|
save(sprintf(strcat(this.SaveDirectory, '/Run_%03i/psi_gs.mat'),Params.njob),'psi','muchem','Observ','Transf','Params','VDk','V');
|
|
disp('Completed and saved!');
|
|
end
|
|
end
|
|
|
|
%% Modules
|
|
|
|
function J = compute_gradient(psi, Params, Transf, VDk, V)
|
|
% === Precompute quantities ===
|
|
absPsi2 = abs(psi).^2; % Density |ψ|^2
|
|
absPsi3 = abs(psi).^3; % |ψ|^3
|
|
N = Params.N;
|
|
|
|
% === Kinetic energy operator: -∇²/2 ===
|
|
KEop = 0.5 * (Transf.KX.^2 + Transf.KY.^2 + Transf.KZ.^2);
|
|
HKin = @(w) ifftn(KEop .* fftn(w));
|
|
|
|
% === External potential (already scaled by ℏω₀) ===
|
|
HV = @(w) V .* w;
|
|
|
|
% === Contact interaction: C|ψ|² ===
|
|
C = Params.gs * N;
|
|
Hint = @(w) (C * absPsi2) .* w;
|
|
|
|
% === Dipole-dipole interaction: D U_dd * |ψ|² ===
|
|
frho = fftn(absPsi2);
|
|
Phi_dd = ifftn(frho .* VDk);
|
|
D = Params.gdd * N;
|
|
Hddi = @(w) (D * Phi_dd) .* w;
|
|
|
|
% === Quantum fluctuations: Q|ψ|³ ===
|
|
Q = Params.gammaQF * N^(3/2);
|
|
Hqf = @(w) (Q * absPsi3) .* w;
|
|
|
|
% === Total Hamiltonian operator ===
|
|
H = @(w) HKin(w) + HV(w) + Hint(w) + Hddi(w) + Hqf(w);
|
|
|
|
% === Compute gradient ===
|
|
J = H(psi);
|
|
end
|
|
|
|
function g = compute_g(psi, p, Params, VDk)
|
|
% === Precompute necessary quantities ===
|
|
rho = real(psi .* conj(p)); % Mixed density
|
|
absPsi = abs(psi); % Magnitude of wavefunction
|
|
N = Params.N;
|
|
|
|
% === Contact interaction term ===
|
|
C = Params.gs * N;
|
|
g_contact = @(w) (C * rho) .* w;
|
|
|
|
% === Dipole-dipole interaction term ===
|
|
rhotilde = fftn(rho);
|
|
Phi_dd = ifftn(rhotilde .* VDk); % Convolution with DDI kernel
|
|
D = Params.gdd * N;
|
|
g_ddi = @(w) (D * Phi_dd) .* w;
|
|
|
|
% === Quantum fluctuation (LHY) correction ===
|
|
Q = Params.gammaQF * N^(3/2); % LHY prefactor
|
|
g_qf = @(w) ((3/2) * Q * absPsi .* rho) .* w;
|
|
|
|
% === Total effective operator applied to psi ===
|
|
g_total = @(w) g_contact(w) + g_ddi(w) + g_qf(w);
|
|
|
|
% === Multiply by factor of 2, as per gradient expression ===
|
|
g = 2 * g_total(psi);
|
|
end
|
|
|
|
% Optimal direction
|
|
function theta = compute_optimal_theta(p, muchem, psi, Params, Transf, VDk, V)
|
|
|
|
Hpsi = compute_gradient(psi, Params, Transf, VDk, V);
|
|
Hp = compute_gradient(p, Params, Transf, VDk, V);
|
|
g = compute_g(psi, p, Params, VDk);
|
|
gamma = 1 / sqrt(inner_product(p, p, Transf));
|
|
numerator = gamma * real(inner_product(p, Hpsi, Transf));
|
|
denominator = muchem - (gamma^2 * (real(inner_product(p, Hp, Transf)) + real(inner_product(g, p, Transf))));
|
|
|
|
theta = numerator / denominator;
|
|
end
|
|
|
|
% Optimal step size via Polak-Ribiere
|
|
function beta = compute_beta(residual_new, residual_old, Transf)
|
|
beta = inner_product(residual_new, (residual_new - residual_old), Transf) / inner_product(residual_old, residual_old, Transf);
|
|
beta = max(0, real(beta));
|
|
end
|
|
|
|
function s = inner_product(u, v, Transf)
|
|
s = sum(conj(u(:)) .* v(:)) * Transf.dx * Transf.dy * Transf.dz;
|
|
end
|