From e0991a9a29b914bf2b30289e1397b1c26aa28e91 Mon Sep 17 00:00:00 2001 From: Karthik Chandrashekara Date: Thu, 3 Apr 2025 17:17:34 +0200 Subject: [PATCH] Latest script - includes major modification to gradient descent function and minor aesthetic changes. --- .../+Plotter/visualizeGSWavefunction.m | 14 +- Dipolar-Gas-Simulator/+Scripts/run_locally.m | 25 +- .../+Scripts/run_on_cluster.m | 56 +++- .../+Simulator/@DipolarGas/DipolarGas.m | 38 ++- .../@DipolarGas/runGradientDescent.m | 294 +++++++++++++----- .../@DipolarGas/setupWavefunction.m | 38 +-- 6 files changed, 324 insertions(+), 141 deletions(-) diff --git a/Dipolar-Gas-Simulator/+Plotter/visualizeGSWavefunction.m b/Dipolar-Gas-Simulator/+Plotter/visualizeGSWavefunction.m index 1bb9ef4..99dd50e 100644 --- a/Dipolar-Gas-Simulator/+Plotter/visualizeGSWavefunction.m +++ b/Dipolar-Gas-Simulator/+Plotter/visualizeGSWavefunction.m @@ -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; diff --git a/Dipolar-Gas-Simulator/+Scripts/run_locally.m b/Dipolar-Gas-Simulator/+Scripts/run_locally.m index 8dc0c14..40e7774 100644 --- a/Dipolar-Gas-Simulator/+Scripts/run_locally.m +++ b/Dipolar-Gas-Simulator/+Scripts/run_locally.m @@ -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) diff --git a/Dipolar-Gas-Simulator/+Scripts/run_on_cluster.m b/Dipolar-Gas-Simulator/+Scripts/run_on_cluster.m index 85e3981..5280b0e 100644 --- a/Dipolar-Gas-Simulator/+Scripts/run_on_cluster.m +++ b/Dipolar-Gas-Simulator/+Scripts/run_on_cluster.m @@ -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(); \ No newline at end of file +[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(); +%} \ No newline at end of file diff --git a/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/DipolarGas.m b/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/DipolarGas.m index 07384c2..288d30d 100644 --- a/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/DipolarGas.m +++ b/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/DipolarGas.m @@ -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; diff --git a/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/runGradientDescent.m b/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/runGradientDescent.m index b9a5301..8ac935c 100644 --- a/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/runGradientDescent.m +++ b/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/runGradientDescent.m @@ -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 \ No newline at end of file diff --git a/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/setupWavefunction.m b/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/setupWavefunction.m index 973e017..2b6dc7c 100644 --- a/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/setupWavefunction.m +++ b/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/setupWavefunction.m @@ -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 \ No newline at end of file