Latest script - includes major modification to gradient descent function and minor aesthetic changes.
This commit is contained in:
parent
767cec6d86
commit
e0991a9a29
@ -11,18 +11,10 @@ function visualizeGSWavefunction(folder_path, run_index)
|
||||
Params = Data.Params;
|
||||
Transf = Data.Transf;
|
||||
Observ = Data.Observ;
|
||||
if isgpuarray(Data.psi)
|
||||
psi = gather(Data.psi);
|
||||
else
|
||||
psi = Data.psi;
|
||||
end
|
||||
|
||||
if isgpuarray(Data.Observ.residual)
|
||||
Observ.residual = gather(Data.Observ.residual);
|
||||
else
|
||||
Observ.residual = Data.Observ.residual;
|
||||
end
|
||||
|
||||
if isgpuarray(Data.psi), psi = gather(Data.psi); end
|
||||
if isgpuarray(Data.Observ.residual), Observ.residual = gather(Data.Observ.residual); end
|
||||
|
||||
% Axes scaling and coordinates in micrometers
|
||||
x = Transf.x * Params.l0 * 1e6;
|
||||
y = Transf.y * Params.l0 * 1e6;
|
||||
|
||||
@ -6,7 +6,7 @@
|
||||
|
||||
OptionsStruct = struct;
|
||||
|
||||
OptionsStruct.NumberOfAtoms = 8E4;
|
||||
OptionsStruct.NumberOfAtoms = 40000;
|
||||
OptionsStruct.DipolarPolarAngle = deg2rad(0);
|
||||
OptionsStruct.DipolarAzimuthAngle = 0;
|
||||
OptionsStruct.ScatteringLength = 95;
|
||||
@ -14,19 +14,20 @@ OptionsStruct.ScatteringLength = 95;
|
||||
OptionsStruct.TrapFrequencies = [30, 60, 90];
|
||||
OptionsStruct.TrapPotentialType = 'Harmonic';
|
||||
|
||||
OptionsStruct.NumberOfGridPoints = [256, 128, 128];
|
||||
OptionsStruct.NumberOfGridPoints = [128, 64, 64];
|
||||
OptionsStruct.Dimensions = [30, 20, 20];
|
||||
OptionsStruct.UseApproximationForLHY = true;
|
||||
OptionsStruct.IncludeDDICutOff = true;
|
||||
OptionsStruct.CutoffType = 'Cylindrical';
|
||||
OptionsStruct.SimulationMode = 'EnergyMinimization'; % 'ImaginaryTimeEvolution' | 'RealTimeEvolution' | 'EnergyMinimization'
|
||||
OptionsStruct.MaxIterationsForGD = 2E4;
|
||||
% OptionsStruct.TimeStepSize = 1E-4; % in s
|
||||
% OptionsStruct.MinimumTimeStepSize = 2E-10; % in s
|
||||
% OptionsStruct.TimeCutOff = 2E6; % in s
|
||||
% OptionsStruct.EnergyTolerance = 5E-10;
|
||||
% OptionsStruct.ResidualTolerance = 1E-08;
|
||||
OptionsStruct.NoiseScaleFactor = 0.01;
|
||||
OptionsStruct.GradientDescentMethod = 'HeavyBall'; % 'HeavyBall' | 'NonLinearCGD'
|
||||
OptionsStruct.MaxIterationsForGD = 2E5;
|
||||
OptionsStruct.TimeStepSize = 1E-4; % in s
|
||||
OptionsStruct.MinimumTimeStepSize = 2E-10; % in s
|
||||
OptionsStruct.TimeCutOff = 2E6; % in s
|
||||
OptionsStruct.EnergyTolerance = 5E-10;
|
||||
OptionsStruct.ResidualTolerance = 1E-08;
|
||||
OptionsStruct.NoiseScaleFactor = 0.010;
|
||||
|
||||
OptionsStruct.PlotLive = true;
|
||||
OptionsStruct.JobNumber = 0;
|
||||
@ -66,7 +67,7 @@ SaveDirectory = './Results/Data_3D/ApproximateLHY/AspectRatio2_8';
|
||||
% SaveDirectory = './Results/Data_3D/ApproximateLHY/AspectRatio3_7';
|
||||
% SaveDirectory = './Results/Data_3D/ApproximateLHY/AspectRatio3_8';
|
||||
% SaveDirectory = './Results/Data_3D/ApproximateLHY/AspectRatio3_9';
|
||||
JobNumber = 2;
|
||||
JobNumber = 3;
|
||||
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)
|
||||
%%
|
||||
% SaveDirectory = './Results/Data_3D/ApproximateLHY/BeyondSSD_SSD';
|
||||
@ -75,6 +76,10 @@ SaveDirectory = './Results/Data_3D/ApproximateLHY/BeyondSSD_Honeycomb';
|
||||
JobNumber = 0;
|
||||
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)
|
||||
%%
|
||||
SaveDirectory = './Results/Data_3D/GradientDescent';
|
||||
JobNumber = 0;
|
||||
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)
|
||||
%%
|
||||
|
||||
% To reproduce results from the Blair Blakie paper:
|
||||
% (n*add^2, as/add)
|
||||
|
||||
@ -1,13 +1,13 @@
|
||||
%% - N = 8E4
|
||||
%% Scaled parameters
|
||||
|
||||
ScalingFactor = (0.8/5)^2;
|
||||
ScalingFactor = (4/5)^2;
|
||||
|
||||
OptionsStruct = struct;
|
||||
|
||||
OptionsStruct.NumberOfAtoms = sqrt(ScalingFactor) * 5E5;
|
||||
OptionsStruct.DipolarPolarAngle = deg2rad(0);
|
||||
OptionsStruct.DipolarAzimuthAngle = 0;
|
||||
OptionsStruct.ScatteringLength = 85;
|
||||
OptionsStruct.ScatteringLength = 75;
|
||||
|
||||
AspectRatio = 2.8;
|
||||
HorizontalTrapFrequency = 125/ScalingFactor;
|
||||
@ -16,7 +16,7 @@ OptionsStruct.TrapFrequencies = [HorizontalTrapFrequency, HorizontalTrap
|
||||
OptionsStruct.TrapPotentialType = 'Harmonic';
|
||||
|
||||
OptionsStruct.NumberOfGridPoints = [128, 128, 64];
|
||||
OptionsStruct.Dimensions = [4, 4, 4];
|
||||
OptionsStruct.Dimensions = [18, 18, 18];
|
||||
OptionsStruct.UseApproximationForLHY = true;
|
||||
OptionsStruct.IncludeDDICutOff = true;
|
||||
OptionsStruct.CutoffType = 'Cylindrical';
|
||||
@ -29,7 +29,7 @@ OptionsStruct.ResidualTolerance = 1E-08;
|
||||
OptionsStruct.NoiseScaleFactor = 0.01;
|
||||
|
||||
OptionsStruct.PlotLive = false;
|
||||
OptionsStruct.JobNumber = 0;
|
||||
OptionsStruct.JobNumber = 2;
|
||||
OptionsStruct.RunOnGPU = true;
|
||||
OptionsStruct.SaveData = true;
|
||||
OptionsStruct.SaveDirectory = sprintf('./Results/Data_3D/ApproximateLHY/AspectRatio%s', strrep(num2str(AspectRatio), '.', '_'));
|
||||
@ -41,4 +41,48 @@ pot = Simulator.Potentials(options{:});
|
||||
sim.Potential = pot.trap();
|
||||
|
||||
%-% Run Simulation %-%
|
||||
[Params, Transf, psi, V, VDk] = sim.run();
|
||||
[Params, Transf, psi, V, VDk] = sim.run();
|
||||
|
||||
|
||||
%{
|
||||
%% - Gradient Descent Test
|
||||
|
||||
OptionsStruct = struct;
|
||||
|
||||
OptionsStruct.NumberOfAtoms = 8E4;
|
||||
OptionsStruct.DipolarPolarAngle = deg2rad(0);
|
||||
OptionsStruct.DipolarAzimuthAngle = 0;
|
||||
OptionsStruct.ScatteringLength = 95;
|
||||
|
||||
OptionsStruct.TrapFrequencies = [30, 60, 90];
|
||||
OptionsStruct.TrapPotentialType = 'Harmonic';
|
||||
|
||||
OptionsStruct.NumberOfGridPoints = [256, 128, 128];
|
||||
OptionsStruct.Dimensions = [30, 20, 20];
|
||||
OptionsStruct.UseApproximationForLHY = true;
|
||||
OptionsStruct.IncludeDDICutOff = true;
|
||||
OptionsStruct.CutoffType = 'Cylindrical';
|
||||
OptionsStruct.SimulationMode = 'EnergyMinimization'; % 'ImaginaryTimeEvolution' | 'RealTimeEvolution' | 'EnergyMinimization'
|
||||
OptionsStruct.MaxIterationsForGD = 2E4;
|
||||
% OptionsStruct.TimeStepSize = 1E-4; % in s
|
||||
% OptionsStruct.MinimumTimeStepSize = 2E-10; % in s
|
||||
% OptionsStruct.TimeCutOff = 2E6; % in s
|
||||
% OptionsStruct.EnergyTolerance = 5E-10;
|
||||
% OptionsStruct.ResidualTolerance = 1E-08;
|
||||
OptionsStruct.NoiseScaleFactor = 0.01;
|
||||
|
||||
OptionsStruct.PlotLive = true;
|
||||
OptionsStruct.JobNumber = 0;
|
||||
OptionsStruct.RunOnGPU = false;
|
||||
OptionsStruct.SaveData = true;
|
||||
OptionsStruct.SaveDirectory = './Results/Data_3D/GradientDescent';
|
||||
options = Helper.convertstruct2cell(OptionsStruct);
|
||||
clear OptionsStruct
|
||||
|
||||
sim = Simulator.DipolarGas(options{:});
|
||||
pot = Simulator.Potentials(options{:});
|
||||
sim.Potential = pot.trap(); % + pot.repulsive_chopstick();
|
||||
|
||||
%-% Run Simulation %-%
|
||||
[Params, Transf, psi, V, VDk] = sim.run();
|
||||
%}
|
||||
@ -17,6 +17,7 @@ classdef DipolarGas < handle & matlab.mixin.Copyable
|
||||
EnergyTolerance;
|
||||
ResidualTolerance;
|
||||
NoiseScaleFactor;
|
||||
GradientDescentMethod;
|
||||
MaxIterationsForGD;
|
||||
|
||||
Calculator;
|
||||
@ -69,7 +70,9 @@ classdef DipolarGas < handle & matlab.mixin.Copyable
|
||||
addParameter(p, 'ResidualTolerance', 1e-10,...
|
||||
@(x) assert(isnumeric(x) && isscalar(x) && (x > 0)));
|
||||
addParameter(p, 'NoiseScaleFactor', 4,...
|
||||
@(x) assert(isnumeric(x) && isscalar(x) && (x > 0)));
|
||||
@(x) assert(isnumeric(x) && isscalar(x) && (x >= 0)));
|
||||
addParameter(p, 'GradientDescentMethod', 'HeavyBall',...
|
||||
@(x) assert(any(strcmpi(x,{'HeavyBall','NonLinearCGD'}))));
|
||||
addParameter(p, 'MaxIterationsForGD', 100,...
|
||||
@(x) assert(isnumeric(x) && isscalar(x) && (x > 0)));
|
||||
|
||||
@ -92,22 +95,23 @@ classdef DipolarGas < handle & matlab.mixin.Copyable
|
||||
|
||||
p.parse(varargin{:});
|
||||
|
||||
this.NumberOfAtoms = p.Results.NumberOfAtoms;
|
||||
this.DipolarPolarAngle = p.Results.DipolarPolarAngle;
|
||||
this.DipolarAzimuthAngle = p.Results.DipolarAzimuthAngle;
|
||||
this.ScatteringLength = p.Results.ScatteringLength;
|
||||
this.TrapFrequencies = p.Results.TrapFrequencies;
|
||||
this.NumberOfGridPoints = p.Results.NumberOfGridPoints;
|
||||
this.Dimensions = p.Results.Dimensions;
|
||||
this.Potential = NaN;
|
||||
this.SimulationMode = p.Results.SimulationMode;
|
||||
this.TimeStepSize = p.Results.TimeStepSize;
|
||||
this.MinimumTimeStepSize = p.Results.MinimumTimeStepSize;
|
||||
this.TimeCutOff = p.Results.TimeCutOff;
|
||||
this.EnergyTolerance = p.Results.EnergyTolerance;
|
||||
this.ResidualTolerance = p.Results.ResidualTolerance;
|
||||
this.NoiseScaleFactor = p.Results.NoiseScaleFactor;
|
||||
this.MaxIterationsForGD = p.Results.MaxIterationsForGD;
|
||||
this.NumberOfAtoms = p.Results.NumberOfAtoms;
|
||||
this.DipolarPolarAngle = p.Results.DipolarPolarAngle;
|
||||
this.DipolarAzimuthAngle = p.Results.DipolarAzimuthAngle;
|
||||
this.ScatteringLength = p.Results.ScatteringLength;
|
||||
this.TrapFrequencies = p.Results.TrapFrequencies;
|
||||
this.NumberOfGridPoints = p.Results.NumberOfGridPoints;
|
||||
this.Dimensions = p.Results.Dimensions;
|
||||
this.Potential = NaN;
|
||||
this.SimulationMode = p.Results.SimulationMode;
|
||||
this.TimeStepSize = p.Results.TimeStepSize;
|
||||
this.MinimumTimeStepSize = p.Results.MinimumTimeStepSize;
|
||||
this.TimeCutOff = p.Results.TimeCutOff;
|
||||
this.EnergyTolerance = p.Results.EnergyTolerance;
|
||||
this.ResidualTolerance = p.Results.ResidualTolerance;
|
||||
this.NoiseScaleFactor = p.Results.NoiseScaleFactor;
|
||||
this.GradientDescentMethod = p.Results.GradientDescentMethod;
|
||||
this.MaxIterationsForGD = p.Results.MaxIterationsForGD;
|
||||
|
||||
this.IncludeDDICutOff = p.Results.IncludeDDICutOff;
|
||||
this.UseApproximationForLHY = p.Results.UseApproximationForLHY;
|
||||
|
||||
@ -1,90 +1,190 @@
|
||||
function [psi] = runGradientDescent(this,psi,Params,Transf,VDk,V,Observ)
|
||||
|
||||
format long;
|
||||
|
||||
% Convergence Criteria:
|
||||
epsilon = 1E-6;
|
||||
alpha = 1E-3;
|
||||
beta = 0.9;
|
||||
|
||||
Observ.residual = 1;
|
||||
Observ.res = 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
|
||||
switch this.GradientDescentMethod
|
||||
case 'HeavyBall'
|
||||
% Convergence Criteria:
|
||||
epsilon = 1E-6;
|
||||
alpha = 1E-3;
|
||||
beta = 0.9;
|
||||
|
||||
% 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;
|
||||
else
|
||||
|
||||
% Update psi using heavy-ball method
|
||||
psi_new = (1 + beta) * psi - alpha * J - beta * psi_old;
|
||||
psi_old = psi;
|
||||
Observ.residual = 1;
|
||||
Observ.res = 1;
|
||||
psi_old = psi; % Previous psi value (for heavy-ball method)
|
||||
|
||||
% Normalize psi
|
||||
Norm = sum(abs(psi_new(:)).^2) * Transf.dx * Transf.dy * Transf.dz;
|
||||
psi = sqrt(Params.N) * 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);
|
||||
E = E/Norm;
|
||||
Observ.EVec = [Observ.EVec E];
|
||||
|
||||
% Collect chemical potentials
|
||||
Observ.mucVec = [Observ.mucVec muchem];
|
||||
|
||||
% Collect residuals
|
||||
Observ.residual = [Observ.residual residual];
|
||||
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');
|
||||
if this.PlotLive
|
||||
Plotter.plotLive(psi,Params,Transf,Observ)
|
||||
drawnow
|
||||
end
|
||||
end
|
||||
|
||||
% Minimization Loop
|
||||
for idx = 1:this.MaxIterationsForGD
|
||||
|
||||
% Compute gradient
|
||||
J = compute_gradient(psi, Params, Transf, VDk, V);
|
||||
|
||||
% Calculate chemical potential and norm
|
||||
% Can also be calculated as --> muchem = this.Calculator.calculateChemicalPotential(psi,Params,Transf,VDk,V);
|
||||
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;
|
||||
else
|
||||
|
||||
% 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);
|
||||
|
||||
% Write output at specified intervals
|
||||
if mod(idx,100) == 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 residual];
|
||||
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
|
||||
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!');
|
||||
case 'NonLinearCGD'
|
||||
% Define the function handle
|
||||
f = @(X) this.Calculator.calculateTotalEnergy(X, Params, Transf, VDk, V)/Params.N;
|
||||
|
||||
% Convergence Criteria:
|
||||
Epsilon = 1E-5;
|
||||
|
||||
% Iteration Counter:
|
||||
i = 1;
|
||||
Observ.residual = 1;
|
||||
Observ.res = 1;
|
||||
|
||||
% Initialize the PrematureExitFlag to false
|
||||
PrematureExitFlag = false;
|
||||
|
||||
if this.PlotLive
|
||||
Plotter.plotLive(psi,Params,Transf,Observ)
|
||||
drawnow
|
||||
end
|
||||
|
||||
% Minimization Loop
|
||||
while true
|
||||
% Compute gradient
|
||||
J = compute_gradient(psi, Params, Transf, VDk, V);
|
||||
|
||||
% Check stopping criterion (Gradient norm)
|
||||
if norm(J(:)) < Epsilon
|
||||
disp('Tolerance reached: Gradient norm is below the specified epsilon.');
|
||||
PrematureExitFlag = true; % Set flag to indicate premature exit
|
||||
break;
|
||||
elseif i >= this.MaxIterationsForCGD
|
||||
disp('Maximum number of iterations for CGD reached.');
|
||||
PrematureExitFlag = true; % Set flag to indicate premature exit
|
||||
break;
|
||||
end
|
||||
|
||||
% Initialize search direction if first iteration
|
||||
if i == 1
|
||||
S = -J;
|
||||
else
|
||||
% Update search direction
|
||||
S = update_search_direction(S, J, J_old);
|
||||
end
|
||||
|
||||
% Step Size Optimization (Line Search)
|
||||
alpha = optimize_step_size(f, psi, S, Params, Transf, VDk, V);
|
||||
|
||||
% Update solution
|
||||
psi = psi + alpha * S;
|
||||
|
||||
% Normalize psi
|
||||
Norm = sum(abs(psi(:)).^2) * Transf.dx * Transf.dy * Transf.dz;
|
||||
psi = sqrt(Params.N) * psi / sqrt(Norm);
|
||||
|
||||
% Store old gradient
|
||||
J_old = J;
|
||||
i = i + 1;
|
||||
muchem = this.Calculator.calculateChemicalPotential(psi,Params,Transf,VDk,V);
|
||||
|
||||
if mod(i,500) == 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)
|
||||
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);
|
||||
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
|
||||
|
||||
% 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
|
||||
%% Modules
|
||||
% Numerical Gradient Calculation using the finite differences method
|
||||
function J = compute_gradient(psi, Params, Transf, VDk, V)
|
||||
|
||||
@ -112,4 +212,42 @@ function J = compute_gradient(psi, Params, Transf, VDk, V)
|
||||
|
||||
J = H(psi);
|
||||
|
||||
end
|
||||
|
||||
% Backtracking Line Search (Step Size Optimization)
|
||||
function alpha = optimize_step_size(f, X, S, Params, Transf, VDk, V)
|
||||
alpha = 1; % Initial step size
|
||||
rho = 0.5; % Step size reduction factor
|
||||
c = 1E-4; % Armijo condition constant
|
||||
max_iter = 100; % Max iterations for backtracking
|
||||
tol = 1E-4; % Tolerance for stopping
|
||||
|
||||
grad = compute_gradient(X, Params, Transf, VDk, V); % Compute gradient once
|
||||
f_X = f(X); % Evaluate f(X) once
|
||||
|
||||
for k = 1:max_iter
|
||||
% Evaluate Armijo condition with precomputed f(X) and grad
|
||||
if f(X + alpha * S) <= f_X + c * alpha * (S(:)' * grad(:))
|
||||
break;
|
||||
else
|
||||
alpha = rho * alpha; % Reduce the step size
|
||||
end
|
||||
|
||||
% Early stopping if step size becomes too small
|
||||
if alpha < tol
|
||||
break;
|
||||
end
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
% Update Search Direction
|
||||
function S_new = update_search_direction(S, J_new, J_old)
|
||||
% (Fletcher-Reeves method)
|
||||
% beta = (norm(J_new(:))^2) / (norm(J_old(:))^2);
|
||||
% S_new = -J_new + beta * S;
|
||||
|
||||
% (Polak-Ribiere method)
|
||||
beta = max(0, (J_new(:)' * (J_new(:) - J_old(:))) / (norm(J_old(:))^2));
|
||||
S_new = -J_new + beta * S;
|
||||
end
|
||||
@ -1,32 +1,32 @@
|
||||
function [psi] = setupWavefunction(~,Params,Transf)
|
||||
|
||||
X = Transf.X;
|
||||
Y = Transf.Y;
|
||||
Z = Transf.Z;
|
||||
|
||||
ellx = sqrt(Params.hbar/(Params.m*Params.wx))/Params.l0;
|
||||
elly = sqrt(Params.hbar/(Params.m*Params.wy))/Params.l0;
|
||||
ellz = sqrt(Params.hbar/(Params.m*Params.wz))/Params.l0;
|
||||
X = Transf.X;
|
||||
Y = Transf.Y;
|
||||
Z = Transf.Z;
|
||||
|
||||
Rx = 8*ellx;
|
||||
Ry = 8*elly;
|
||||
Rz = 8*ellz;
|
||||
X0 = 0.0*Transf.Xmax;
|
||||
Y0 = 0.0*Transf.Ymax;
|
||||
Z0 = 0*Transf.Zmax;
|
||||
ellx = sqrt(Params.hbar/(Params.m*Params.wx))/Params.l0;
|
||||
elly = sqrt(Params.hbar/(Params.m*Params.wy))/Params.l0;
|
||||
ellz = sqrt(Params.hbar/(Params.m*Params.wz))/Params.l0;
|
||||
|
||||
Rx = 4.0*ellx;
|
||||
Ry = 4.0*elly;
|
||||
Rz = 4.0*ellz;
|
||||
X0 = 0.0*Transf.Xmax;
|
||||
Y0 = 0.0*Transf.Ymax;
|
||||
Z0 = 0.0*Transf.Zmax;
|
||||
|
||||
psi = exp(-(X-X0).^2/Rx^2-(Y-Y0).^2/Ry^2-(Z-Z0).^2/Rz^2);
|
||||
cur_norm = sum(abs(psi(:)).^2)*Transf.dx*Transf.dy*Transf.dz;
|
||||
psi = psi/sqrt(cur_norm);
|
||||
|
||||
% --- Adding some noise ---
|
||||
r = normrnd(0,1,size(X));
|
||||
theta = rand(size(X));
|
||||
noise = r.*exp(2*pi*1i*theta);
|
||||
psi = psi + Params.nsf*noise;
|
||||
r = normrnd(0,1,size(X));
|
||||
theta = rand(size(X));
|
||||
noise = r.*exp(2*pi*1i*theta);
|
||||
psi = psi + Params.nsf*noise;
|
||||
|
||||
% Renormalize wavefunction
|
||||
Norm = sum(abs(psi(:)).^2)*Transf.dx*Transf.dy*Transf.dz;
|
||||
psi = sqrt(Params.N)*psi/sqrt(Norm);
|
||||
Norm = sum(abs(psi(:)).^2)*Transf.dx*Transf.dy*Transf.dz;
|
||||
psi = sqrt(Params.N)*psi/sqrt(Norm);
|
||||
|
||||
end
|
||||
Loading…
Reference in New Issue
Block a user