karthik/Calculations
1
0
Fork 0
Code Issues Pull Requests repo.project_board Releases Wiki Activity

Added functionality to load a initial wavefunction if necessary to bias the solver, modified the conjugate gradient descent algorithm.

Browse Source
This commit is contained in:
Karthik 2025-04-04 21:26:25 +02:00
parent e0991a9a29
commit 6566f53559
4 changed files with 153 additions and 158 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

69
Dipolar-Gas-Simulator/+Scripts/run_locally.m
View File

@ -20,7 +20,7 @@ OptionsStruct.UseApproximationForLHY = true;
OptionsStruct.IncludeDDICutOff = true; OptionsStruct.IncludeDDICutOff = true;
OptionsStruct.CutoffType = 'Cylindrical'; OptionsStruct.CutoffType = 'Cylindrical';
OptionsStruct.SimulationMode = 'EnergyMinimization'; % 'ImaginaryTimeEvolution' | 'RealTimeEvolution' | 'EnergyMinimization' OptionsStruct.SimulationMode = 'EnergyMinimization'; % 'ImaginaryTimeEvolution' | 'RealTimeEvolution' | 'EnergyMinimization'
OptionsStruct.GradientDescentMethod = 'HeavyBall'; % 'HeavyBall' | 'NonLinearCGD' OptionsStruct.GradientDescentMethod = 'NonLinearCGD'; % 'HeavyBall' | 'NonLinearCGD'
OptionsStruct.MaxIterationsForGD = 2E5; OptionsStruct.MaxIterationsForGD = 2E5;
OptionsStruct.TimeStepSize = 1E-4; % in s OptionsStruct.TimeStepSize = 1E-4; % in s
OptionsStruct.MinimumTimeStepSize = 2E-10; % in s OptionsStruct.MinimumTimeStepSize = 2E-10; % in s
@ -50,35 +50,6 @@ sim.Potential = pot.trap(); % + pot.repulsive_chopstick(
Plotter.visualizeTrapPotential(sim.Potential,Params,Transf) Plotter.visualizeTrapPotential(sim.Potential,Params,Transf)
%% - Plot initial wavefunction %% - Plot initial wavefunction
Plotter.visualizeWavefunction(psi,Params,Transf) Plotter.visualizeWavefunction(psi,Params,Transf)
%% - Plot GS wavefunction
% SaveDirectory = './Results/Data_3D/CompleteLHY/AspectRatio2_8';
% SaveDirectory = './Results/Data_3D/CompleteLHY/AspectRatio3_7';
SaveDirectory = './Results/Data_3D/CompleteLHY/AspectRatio3_8';
JobNumber = 0;
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)
%%
% SaveDirectory = './Results/Data_3D/CompleteLHY/BeyondSSD_SSD';
% SaveDirectory = './Results/Data_3D/CompleteLHY/BeyondSSD_Labyrinth';
SaveDirectory = './Results/Data_3D/CompleteLHY/BeyondSSD_Honeycomb';
JobNumber = 0;
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)
%% - Plot GS wavefunction
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 = 3;
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)
%%
% SaveDirectory = './Results/Data_3D/ApproximateLHY/BeyondSSD_SSD';
% SaveDirectory = './Results/Data_3D/ApproximateLHY/BeyondSSD_Labyrinth';
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: % To reproduce results from the Blair Blakie paper:
@ -438,4 +409,40 @@ JobNumber = 2; % 82
% JobNumber = 3; % 83 % JobNumber = 3; % 83
Plotter.visualizeGSWavefunction2D(SaveDirectory, JobNumber) Plotter.visualizeGSWavefunction2D(SaveDirectory, JobNumber)
[contrast, period_X, period_Y] = Scripts.analyzeGSWavefunction_in_plane_trap(SaveDirectory, JobNumber); [contrast, period_X, period_Y] = Scripts.analyzeGSWavefunction_in_plane_trap(SaveDirectory, JobNumber);
%% - Plot GS wavefunction
% SaveDirectory = './Results/Data_3D/CompleteLHY/AspectRatio2_8';
% SaveDirectory = './Results/Data_3D/CompleteLHY/AspectRatio3_7';
SaveDirectory = './Results/Data_3D/CompleteLHY/AspectRatio3_8';
JobNumber = 0;
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)
%%
% SaveDirectory = './Results/Data_3D/CompleteLHY/BeyondSSD_SSD';
% SaveDirectory = './Results/Data_3D/CompleteLHY/BeyondSSD_Labyrinth';
SaveDirectory = './Results/Data_3D/CompleteLHY/BeyondSSD_Honeycomb';
JobNumber = 0;
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)
%% - Plot GS wavefunction
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 = 3;
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)
%%
% SaveDirectory = './Results/Data_3D/ApproximateLHY/BeyondSSD_SSD';
% SaveDirectory = './Results/Data_3D/ApproximateLHY/BeyondSSD_Labyrinth';
SaveDirectory = './Results/Data_3D/ApproximateLHY/BeyondSSD_Honeycomb';
JobNumber = 0;
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)
%%
SaveDirectory = './Results/Data_3D/TiltedDipoles15';
JobNumber = 0;
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)
%%
SaveDirectory = './Results/Data_3D/TiltedDipoles30';
JobNumber = 0;
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)
%%
SaveDirectory = './Results/Data_3D/TiltedDipoles45';
JobNumber = 1;
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)

77
Dipolar-Gas-Simulator/+Scripts/run_on_cluster.m
View File

@ -1,18 +1,20 @@
%% Scaled parameters theta = 0;
phi = 0;
ScalingFactor = (4/5)^2; % - SSD: N = 1E5, as = 86ao
OptionsStruct = struct; OptionsStruct = struct;
OptionsStruct.NumberOfAtoms = sqrt(ScalingFactor) * 5E5; OptionsStruct.NumberOfAtoms = 1E5;
OptionsStruct.DipolarPolarAngle = deg2rad(0); OptionsStruct.DipolarPolarAngle = deg2rad(theta);
OptionsStruct.DipolarAzimuthAngle = 0; OptionsStruct.DipolarAzimuthAngle = deg2rad(phi);
OptionsStruct.ScatteringLength = 75; OptionsStruct.ScatteringLength = 86;
AspectRatio = 2.8; % AspectRatio = 2.0;
HorizontalTrapFrequency = 125/ScalingFactor; % HorizontalTrapFrequency = 125;
VerticalTrapFrequency = AspectRatio * HorizontalTrapFrequency; % VerticalTrapFrequency = AspectRatio * HorizontalTrapFrequency;
OptionsStruct.TrapFrequencies = [HorizontalTrapFrequency, HorizontalTrapFrequency, VerticalTrapFrequency]; % OptionsStruct.TrapFrequencies = [HorizontalTrapFrequency, HorizontalTrapFrequency, VerticalTrapFrequency];
OptionsStruct.TrapFrequencies = [150, 150, 300];
OptionsStruct.TrapPotentialType = 'Harmonic'; OptionsStruct.TrapPotentialType = 'Harmonic';
OptionsStruct.NumberOfGridPoints = [128, 128, 64]; OptionsStruct.NumberOfGridPoints = [128, 128, 64];
@ -21,18 +23,18 @@ OptionsStruct.UseApproximationForLHY = true;
OptionsStruct.IncludeDDICutOff = true; OptionsStruct.IncludeDDICutOff = true;
OptionsStruct.CutoffType = 'Cylindrical'; OptionsStruct.CutoffType = 'Cylindrical';
OptionsStruct.SimulationMode = 'ImaginaryTimeEvolution'; % 'ImaginaryTimeEvolution' | 'RealTimeEvolution' OptionsStruct.SimulationMode = 'ImaginaryTimeEvolution'; % 'ImaginaryTimeEvolution' | 'RealTimeEvolution'
OptionsStruct.TimeStepSize = 1E-4; % in s OptionsStruct.TimeStepSize = 1E-3; % in s
OptionsStruct.MinimumTimeStepSize = 2E-6; % in s OptionsStruct.MinimumTimeStepSize = 2E-6; % in s
OptionsStruct.TimeCutOff = 2E6; % in s OptionsStruct.TimeCutOff = 2E6; % in s
OptionsStruct.EnergyTolerance = 5E-10; OptionsStruct.EnergyTolerance = 5E-10;
OptionsStruct.ResidualTolerance = 1E-08; OptionsStruct.ResidualTolerance = 1E-08;
OptionsStruct.NoiseScaleFactor = 0.01; OptionsStruct.NoiseScaleFactor = 0.01;
OptionsStruct.PlotLive = false; OptionsStruct.PlotLive = true;
OptionsStruct.JobNumber = 2; OptionsStruct.JobNumber = 0;
OptionsStruct.RunOnGPU = true; OptionsStruct.RunOnGPU = false;
OptionsStruct.SaveData = true; OptionsStruct.SaveData = true;
OptionsStruct.SaveDirectory = sprintf('./Results/Data_3D/ApproximateLHY/AspectRatio%s', strrep(num2str(AspectRatio), '.', '_')); OptionsStruct.SaveDirectory = './Results/Data_3D/TiltedDipoles0';
options = Helper.convertstruct2cell(OptionsStruct); options = Helper.convertstruct2cell(OptionsStruct);
clear OptionsStruct clear OptionsStruct
@ -43,46 +45,3 @@ sim.Potential = pot.trap();
%-% Run Simulation %-% %-% 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();
%}

13
Dipolar-Gas-Simulator/+Simulator/@DipolarGas/initialize.m
View File

@ -41,8 +41,19 @@ function [psi,V,VDk] = initialize(this,Params,Transf,TransfRad)
end end
% == Setting up the initial wavefunction == % % == Setting up the initial wavefunction == %
psi = this.setupWavefunction(Params,Transf); WavefunctionFile = fullfile(this.SaveDirectory, 'psi_init.mat');
if isfile(WavefunctionFile)
loadedData = load(WavefunctionFile);
if isgpuarray(loadedData.psi)
psi = gather(loadedData.psi);
else
psi = loadedData.psi;
end
else
psi = this.setupWavefunction(Params,Transf);
end
if this.RunOnGPU if this.RunOnGPU
psi = gpuArray(psi); psi = gpuArray(psi);
end end

152
Dipolar-Gas-Simulator/+Simulator/@DipolarGas/runGradientDescent.m
View File

@ -12,6 +12,7 @@ function [psi] = runGradientDescent(this,psi,Params,Transf,VDk,V,Observ)
Observ.res = 1; Observ.res = 1;
psi_old = psi; % Previous psi value (for heavy-ball method) psi_old = psi; % Previous psi value (for heavy-ball method)
% Live Plotter
if this.PlotLive if this.PlotLive
Plotter.plotLive(psi,Params,Transf,Observ) Plotter.plotLive(psi,Params,Transf,Observ)
drawnow drawnow
@ -24,7 +25,6 @@ function [psi] = runGradientDescent(this,psi,Params,Transf,VDk,V,Observ)
J = compute_gradient(psi, Params, Transf, VDk, V); J = compute_gradient(psi, Params, Transf, VDk, V);
% Calculate chemical potential and norm % 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); muchem = sum(real(conj(psi(:)) .* J(:))) / sum(abs(psi(:)).^2);
% Calculate residual and check convergence % Calculate residual and check convergence
@ -83,12 +83,11 @@ function [psi] = runGradientDescent(this,psi,Params,Transf,VDk,V,Observ)
disp('Saving data...'); disp('Saving data...');
save(sprintf(strcat(this.SaveDirectory, '/Run_%03i/psi_gs.mat'),Params.njob),'psi','muchem','Observ','Transf','Params','VDk','V'); save(sprintf(strcat(this.SaveDirectory, '/Run_%03i/psi_gs.mat'),Params.njob),'psi','muchem','Observ','Transf','Params','VDk','V');
disp('Save complete!'); disp('Save complete!');
case 'NonLinearCGD'
% Define the function handle
f = @(X) this.Calculator.calculateTotalEnergy(X, Params, Transf, VDk, V)/Params.N;
case 'NonLinearCGD'
% Convergence Criteria: % Convergence Criteria:
Epsilon = 1E-5; epsilon = 1E-13;
% Iteration Counter: % Iteration Counter:
i = 1; i = 1;
@ -98,65 +97,80 @@ function [psi] = runGradientDescent(this,psi,Params,Transf,VDk,V,Observ)
% Initialize the PrematureExitFlag to false % Initialize the PrematureExitFlag to false
PrematureExitFlag = false; PrematureExitFlag = false;
% Live plotter
if this.PlotLive if this.PlotLive
Plotter.plotLive(psi,Params,Transf,Observ) Plotter.plotLive(psi,Params,Transf,Observ)
drawnow drawnow
end end
% Minimization Loop % Minimization Loop
while true while true
% Compute gradient % Compute gradient
J = compute_gradient(psi, Params, Transf, VDk, V); J = compute_gradient(psi, Params, Transf, VDk, V);
% Calculate chemical potential
% Check stopping criterion (Gradient norm) muchem = real(conj(psi(:))' * J(:)) / norm(psi(:))^2;
if norm(J(:)) < Epsilon % Calculate residual
disp('Tolerance reached: Gradient norm is below the specified epsilon.'); residual = J - (muchem * psi);
% Compute energy difference between the last two saved energy values
if i == 1
energydifference = NaN;
elseif mod(i,100) == 0 && length(Observ.EVec) > 1
energydifference = abs(Observ.EVec(end) - Observ.EVec(end-1));
end
% Convergence check - if energy difference is below set tolerance, then exit
if energydifference <= epsilon
disp('Tolerance reached: Energy difference is below the specified epsilon.');
PrematureExitFlag = true; % Set flag to indicate premature exit PrematureExitFlag = true; % Set flag to indicate premature exit
break; break;
elseif i >= this.MaxIterationsForCGD elseif i >= this.MaxIterationsForGD % If set maximum number of iterations reached, then exit
disp('Maximum number of iterations for CGD reached.'); disp('Maximum number of iterations for CGD reached.');
PrematureExitFlag = true; % Set flag to indicate premature exit PrematureExitFlag = true; % Set flag to indicate premature exit
break; break;
end end
% Initialize search direction if first iteration % Initialize search direction if first iteration
if i == 1 if i == 1
S = -J; d = -residual;
else else % Compute beta via Polak-Ribiere and create new direction
% Update search direction residual_new = residual;
S = update_search_direction(S, J, J_old); beta = compute_beta(residual_new, residual_old);
d = -residual_new + beta * p_old;
end end
% Step Size Optimization (Line Search) residual_old = residual;
alpha = optimize_step_size(f, psi, S, Params, Transf, VDk, V); p = d - (real(conj(d(:))' * psi(:)) .* 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 % Update solution
psi = psi + alpha * S; psi = (cos(theta).*psi) + (sin(theta).*(p / norm(p(:))));
% Normalize psi % Normalize psi
Norm = sum(abs(psi(:)).^2) * Transf.dx * Transf.dy * Transf.dz; Norm = sum(abs(psi(:)).^2) * Transf.dx * Transf.dy * Transf.dz;
psi = sqrt(Params.N) * psi / sqrt(Norm); psi = sqrt(Params.N) * psi / sqrt(Norm);
i = i + 1;
% Store old gradient
J_old = J; % Calculate chemical potential with new psi
i = i + 1; muchem = real(conj(psi(:))' * J(:)) / norm(psi(:))^2;
muchem = this.Calculator.calculateChemicalPotential(psi,Params,Transf,VDk,V);
if mod(i,100) == 0
if mod(i,500) == 0
% Change in Energy % Collect Energy value
E = this.Calculator.calculateTotalEnergy(psi,Params,Transf,VDk,V); E = this.Calculator.calculateTotalEnergy(psi,Params,Transf,VDk,V);
E = E/Norm; E = E/Norm;
Observ.EVec = [Observ.EVec E]; Observ.EVec = [Observ.EVec E];
% Chemical potential % Collect Chemical potential value
Observ.mucVec = [Observ.mucVec muchem]; Observ.mucVec = [Observ.mucVec muchem];
% Normalized residuals % Collect Normalized residuals
res = this.Calculator.calculateNormalizedResiduals(psi,Params,Transf,VDk,V,muchem); res = this.Calculator.calculateNormalizedResiduals(psi,Params,Transf,VDk,V,muchem);
Observ.residual = [Observ.residual res]; Observ.residual = [Observ.residual res];
Observ.res_idx = Observ.res_idx + 1; Observ.res_idx = Observ.res_idx + 1;
% Live plotter
if this.PlotLive if this.PlotLive
Plotter.plotLive(psi,Params,Transf,Observ) Plotter.plotLive(psi,Params,Transf,Observ)
drawnow drawnow
@ -185,7 +199,7 @@ function [psi] = runGradientDescent(this,psi,Params,Transf,VDk,V,Observ)
end end
%% Modules %% Modules
% Numerical Gradient Calculation using the finite differences method
function J = compute_gradient(psi, Params, Transf, VDk, V) function J = compute_gradient(psi, Params, Transf, VDk, V)
% Operators % Operators
@ -210,44 +224,48 @@ function J = compute_gradient(psi, Params, Transf, VDk, V)
H = @(w) HKin(w) + HV(w) + Hint(w) + Hddi(w) + Hqf(w); H = @(w) HKin(w) + HV(w) + Hint(w) + Hddi(w) + Hqf(w);
J = H(psi); J = H(psi);
end end
% Backtracking Line Search (Step Size Optimization) function g = compute_g(psi, p, Params, VDk)
function alpha = optimize_step_size(f, X, S, Params, Transf, VDk, V)
alpha = 1; % Initial step size rho = real(psi.*conj(p));
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 % Contact interactions
f_X = f(X); % Evaluate f(X) once C = Params.gs*Params.N;
gint = @(w)(C.*rho).*w;
% DDIs
D = Params.gdd*Params.N;
rhotilde = fftn(rho);
Phi = real(ifftn(rhotilde.*VDk));
gaddi = @(w)(D.*Phi).*w;
for k = 1:max_iter % Quantum fluctuations
% Evaluate Armijo condition with precomputed f(X) and grad eps_dd = Params.add/Params.as;
if f(X + alpha * S) <= f_X + c * alpha * (S(:)' * grad(:)) Q = (4/(3*pi^2)) * (C^(5/2)/Params.N) * (1 + ((3*eps_dd^2)/2));
break; gqf = @(w)(((3/2)*Q).*abs(psi).*rho).*w;
else
alpha = rho * alpha; % Reduce the step size gop = @(w) gint(w) + gaddi(w) + gqf(w);
end
g = gop(psi);
% Early stopping if step size becomes too small
if alpha < tol
break;
end
end
end end
% Update Search Direction % Optimal direction via line search
function S_new = update_search_direction(S, J_new, J_old) function theta = compute_optimal_theta(p, muchem, psi, Params, Transf, VDk, V)
% (Fletcher-Reeves method)
% beta = (norm(J_new(:))^2) / (norm(J_old(:))^2); Hpsi = compute_gradient(psi, Params, Transf, VDk, V);
% S_new = -J_new + beta * S; Hp = compute_gradient(p, Params, Transf, VDk, V);
g = compute_g(psi, p, Params, VDk);
numerator = real(conj(p(:))' * Hpsi(:))/norm(p(:));
denominator = muchem - ((conj(p(:))' * Hp(:)) + real(conj(g(:))' * p(:)))/norm(p(:))^2;
% (Polak-Ribiere method) theta = numerator / denominator;
beta = max(0, (J_new(:)' * (J_new(:) - J_old(:))) / (norm(J_old(:))^2)); end
S_new = -J_new + beta * S;
% Optimal step size via Polak-Ribiere
function beta = compute_beta(residual_new, residual_old)
beta = max(0, (residual_new(:)' * (residual_new(:) - residual_old(:))) / (norm(residual_old(:))^2));
end end