From 488c6c9f0bb11308b9e4cc2db2d1e4811210f980 Mon Sep 17 00:00:00 2001 From: Karthik Chandrashekara Date: Wed, 26 Mar 2025 17:20:28 +0100 Subject: [PATCH] Latest updated version of the full 3D simulator. --- Dipolar-Gas-Simulator/+Plotter/plotLive.m | 74 +++-- Dipolar-Gas-Simulator/+Scripts/run_locally.m | 23 +- .../+Scripts/run_on_cluster.m | 264 +++--------------- .../@Calculator/calculateChemicalPotential.m | 36 +-- .../@Calculator/calculateEnergyComponents.m | 35 +-- .../calculateNormalizedResiduals.m | 21 +- .../calculateNumericalHankelTransform.m | 22 +- .../@Calculator/calculateOrderParameter.m | 72 ++--- .../@Calculator/calculatePhaseCoherence.m | 20 +- .../@Calculator/calculateTotalEnergy.m | 38 ++- .../+Simulator/@Calculator/calculateVDk.m | 73 +++++ .../@Calculator/calculateVDkWithCutoff.m | 64 ----- .../+Simulator/@DipolarGas/DipolarGas.m | 36 ++- .../+Simulator/@DipolarGas/initialize.m | 6 +- .../@DipolarGas/propagateWavefunction.m | 108 +++---- .../+Simulator/@DipolarGas/run.m | 15 +- .../+Simulator/@DipolarGas/setupParameters.m | 19 +- .../+Simulator/@DipolarGas/setupSpace.m | 65 +++-- .../@DipolarGas/setupWavefunction.m | 23 +- 19 files changed, 462 insertions(+), 552 deletions(-) create mode 100644 Dipolar-Gas-Simulator/+Simulator/@Calculator/calculateVDk.m delete mode 100644 Dipolar-Gas-Simulator/+Simulator/@Calculator/calculateVDkWithCutoff.m diff --git a/Dipolar-Gas-Simulator/+Plotter/plotLive.m b/Dipolar-Gas-Simulator/+Plotter/plotLive.m index f48288b..cd033be 100644 --- a/Dipolar-Gas-Simulator/+Plotter/plotLive.m +++ b/Dipolar-Gas-Simulator/+Plotter/plotLive.m @@ -3,42 +3,74 @@ function plotLive(psi,Params,Transf,Observ) set(groot, 'defaultAxesTickLabelInterpreter','latex'); set(groot, 'defaultLegendInterpreter','latex'); format long - x = Transf.x*Params.l0*1e6; - y = Transf.y*Params.l0*1e6; - z = Transf.z*Params.l0*1e6; + + % Axes scaling and coordinates in micrometers + x = Transf.x * Params.l0 * 1e6; + y = Transf.y * Params.l0 * 1e6; + z = Transf.z * Params.l0 * 1e6; + dx = x(2)-x(1); dy = y(2)-y(1); dz = z(2)-z(1); - - %Plotting - - subplot(2,3,1) + + % Compute probability density |psi|^2 n = abs(psi).^2; + + %Plotting + figure(1); + clf + set(gcf,'Position', [100, 100, 1600, 900]) + t = tiledlayout(2, 3, 'TileSpacing', 'compact', 'Padding', 'compact'); % 2x3 grid + + nexttile; nxz = squeeze(trapz(n*dy,2)); nyz = squeeze(trapz(n*dx,1)); nxy = squeeze(trapz(n*dz,3)); - plotxz = pcolor(x,z,nxz'); set(plotxz, 'EdgeColor', 'none'); - xlabel('$x$ [$\mu$m]'); ylabel('$z$ [$\mu$m]'); - - subplot(2,3,2) + cbar1 = colorbar; + cbar1.Label.Interpreter = 'latex'; + % cbar1.Ticks = []; % Disable the ticks + colormap(gca, Helper.Colormaps.plasma()) + xlabel('$x$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 14) + ylabel('$z$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 14) + title('$|\Psi(x,y)|^2$', 'Interpreter', 'latex', 'FontSize', 14) + + nexttile; plotyz = pcolor(y,z,nyz'); set(plotyz, 'EdgeColor', 'none'); - xlabel('$y$ [$\mu$m]'); ylabel('$z$ [$\mu$m]'); + cbar1 = colorbar; + cbar1.Label.Interpreter = 'latex'; + % cbar1.Ticks = []; % Disable the ticks + colormap(gca, Helper.Colormaps.plasma()) + xlabel('$y$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 14) + ylabel('$z$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 14) + title('$|\Psi(x,y)|^2$', 'Interpreter', 'latex', 'FontSize', 14) - subplot(2,3,3) + nexttile; plotxy = pcolor(x,y,nxy'); set(plotxy, 'EdgeColor', 'none'); - xlabel('$x$ [$\mu$m]'); ylabel('$y$ [$\mu$m]'); + cbar1 = colorbar; + cbar1.Label.Interpreter = 'latex'; + % cbar1.Ticks = []; % Disable the ticks + colormap(gca, Helper.Colormaps.plasma()) + xlabel('$x$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 14) + ylabel('$y$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 14) + title('$|\Psi(x,y)|^2$', 'Interpreter', 'latex', 'FontSize', 14) - subplot(2,3,4) + nexttile; plot(-log10(Observ.residual),'-b') - ylabel('$-\mathrm{log}_{10}(r)$'); xlabel('steps'); - - subplot(2,3,5) + ylabel('$-\mathrm{log}_{10}(r)$', 'FontSize', 14); xlabel('Time steps', 'FontSize', 14); + title('Residual', 'FontSize', 14); + grid on + + nexttile; plot(Observ.EVec,'-b') - ylabel('$E$'); xlabel('steps'); + ylabel('$E_{tot}$', 'FontSize', 14); xlabel('Time steps', 'FontSize', 14); + title('Total Energy', 'FontSize', 14); + grid on - subplot(2,3,6) + nexttile; plot(Observ.mucVec,'-b') - ylabel('$\mu$'); xlabel('steps'); + ylabel('$\mu$', 'FontSize', 14); xlabel('Time steps', 'FontSize', 14); + title('Chemical Potential', 'FontSize', 14); + grid on end \ No newline at end of file diff --git a/Dipolar-Gas-Simulator/+Scripts/run_locally.m b/Dipolar-Gas-Simulator/+Scripts/run_locally.m index c490ca7..7736548 100644 --- a/Dipolar-Gas-Simulator/+Scripts/run_locally.m +++ b/Dipolar-Gas-Simulator/+Scripts/run_locally.m @@ -6,29 +6,30 @@ OptionsStruct = struct; -OptionsStruct.NumberOfAtoms = 1E5; -OptionsStruct.DipolarPolarAngle = 0; +OptionsStruct.NumberOfAtoms = 5E5; +OptionsStruct.DipolarPolarAngle = deg2rad(0); OptionsStruct.DipolarAzimuthAngle = 0; -OptionsStruct.ScatteringLength = 86; +OptionsStruct.ScatteringLength = 85; -OptionsStruct.TrapFrequencies = [10, 10, 72.4]; -OptionsStruct.TrapDepth = 5; -OptionsStruct.BoxSize = 15; +OptionsStruct.TrapFrequencies = [125, 125, 350]; OptionsStruct.TrapPotentialType = 'Harmonic'; -OptionsStruct.NumberOfGridPoints = [256, 512, 256]; -OptionsStruct.Dimensions = [50, 120, 150]; +OptionsStruct.NumberOfGridPoints = [128, 128, 64]; +OptionsStruct.Dimensions = [30, 30, 15]; OptionsStruct.CutoffType = 'Cylindrical'; OptionsStruct.SimulationMode = 'ImaginaryTimeEvolution'; % 'ImaginaryTimeEvolution' | 'RealTimeEvolution' -OptionsStruct.TimeStepSize = 0.005; % in s -OptionsStruct.NumberOfTimeSteps = 2E6; % in s +OptionsStruct.TimeStepSize = 0.002; % in s +OptionsStruct.MinimumTimeStepSize = 1E-6; % in s +OptionsStruct.TimeCutOff = 1E6; % in s OptionsStruct.EnergyTolerance = 5E-10; OptionsStruct.ResidualTolerance = 1E-05; +OptionsStruct.NoiseScaleFactor = 0.05; +OptionsStruct.PlotLive = true; OptionsStruct.JobNumber = 1; OptionsStruct.RunOnGPU = false; OptionsStruct.SaveData = true; -OptionsStruct.SaveDirectory = './Data'; +OptionsStruct.SaveDirectory = './Results/Data_3D'; options = Helper.convertstruct2cell(OptionsStruct); clear OptionsStruct diff --git a/Dipolar-Gas-Simulator/+Scripts/run_on_cluster.m b/Dipolar-Gas-Simulator/+Scripts/run_on_cluster.m index 1cde2e6..ffa5568 100644 --- a/Dipolar-Gas-Simulator/+Scripts/run_on_cluster.m +++ b/Dipolar-Gas-Simulator/+Scripts/run_on_cluster.m @@ -1,254 +1,82 @@ -%% Tilting of the dipoles -% Atom Number = 1250 ppum -% System size = [sf * unitcell_x, sf * unitcell_x] +%% - Aspect Ratio: 2.8 -ppum = 1250; % Atom Number Density in per micrometers +OptionsStruct = struct; -%% v_z = 500, theta = 0: a_s = 75.00 - -a = 1.8058; -scalingfactor = 5; -lx = scalingfactor*a; -ly = scalingfactor*sqrt(3)*a; - -% Initialize OptionsStruct -OptionsStruct = struct; - -% Assign values to OptionsStruct -OptionsStruct.NumberOfAtoms = ppum * (lx*ly); +OptionsStruct.NumberOfAtoms = 5E5; OptionsStruct.DipolarPolarAngle = deg2rad(0); OptionsStruct.DipolarAzimuthAngle = 0; -OptionsStruct.ScatteringLength = 75.00; +OptionsStruct.ScatteringLength = 85; -OptionsStruct.TrapFrequencies = [0, 0, 500]; -OptionsStruct.TrapPotentialType = 'None'; +AspectRatio = 2.8; +HorizontalTrapFrequency = 125; +VerticalTrapFrequency = AspectRatio * HorizontalTrapFrequency; +OptionsStruct.TrapFrequencies = [HorizontalTrapFrequency, HorizontalTrapFrequency, VerticalTrapFrequency]; +OptionsStruct.TrapPotentialType = 'Harmonic'; -OptionsStruct.NumberOfGridPoints = [128, 128]; -OptionsStruct.Dimensions = [lx, ly]; -OptionsStruct.TimeStepSize = 1E-3; % in s -OptionsStruct.MinimumTimeStepSize = 1E-5; % in s -OptionsStruct.TimeCutOff = 2E6; % in s +OptionsStruct.NumberOfGridPoints = [256, 256, 128]; +OptionsStruct.Dimensions = [30, 30, 15]; +OptionsStruct.CutoffType = 'Cylindrical'; +OptionsStruct.SimulationMode = 'ImaginaryTimeEvolution'; % 'ImaginaryTimeEvolution' | 'RealTimeEvolution' +OptionsStruct.TimeStepSize = 0.002; % in s +OptionsStruct.MinimumTimeStepSize = 1E-6; % in s +OptionsStruct.TimeCutOff = 1E6; % in s OptionsStruct.EnergyTolerance = 5E-10; OptionsStruct.ResidualTolerance = 1E-05; OptionsStruct.NoiseScaleFactor = 0.05; -OptionsStruct.IncludeDDICutOff = false; - -OptionsStruct.MaxIterations = 10; -OptionsStruct.VariationalWidth = 1.15; -OptionsStruct.WidthLowerBound = 0.01; -OptionsStruct.WidthUpperBound = 12; -OptionsStruct.WidthCutoff = 1e-2; OptionsStruct.PlotLive = false; OptionsStruct.JobNumber = 0; OptionsStruct.RunOnGPU = true; OptionsStruct.SaveData = true; -OptionsStruct.SaveDirectory = './Results/Data_TiltingOfDipoles/AdjustedSystemSize/Hz500'; +OptionsStruct.SaveDirectory = sprintf('./Results/Data_3D/AspectRatio%s', strrep(num2str(AspectRatio), '.', '_')); options = Helper.convertstruct2cell(OptionsStruct); clear OptionsStruct -solver = VariationalSolver2D.DipolarGas(options{:}); -pot = VariationalSolver2D.Potentials(options{:}); -solver.Potential = pot.trap(); +sim = Simulator.DipolarGas(options{:}); +pot = Simulator.Potentials(options{:}); +sim.Potential = pot.trap(); % + pot.repulsive_chopstick(); -%-% Run Solver %-% -[Params, Transf, psi, V, VDk] = solver.run(); +%-% Run Simulation %-% +[Params, Transf, psi, V, VDk] = sim.run(); -%% v_z = 500, theta = 5: a_s = 75.00 +%% - Aspect Ratio: 3.7 -a = 1.795; -scalingfactor = 5; -lx = scalingfactor*a; -ly = scalingfactor*sqrt(3)*a; +OptionsStruct = struct; -% Initialize OptionsStruct -OptionsStruct = struct; - -% Assign values to OptionsStruct -OptionsStruct.NumberOfAtoms = ppum * (lx*ly); -OptionsStruct.DipolarPolarAngle = deg2rad(5); +OptionsStruct.NumberOfAtoms = 5E5; +OptionsStruct.DipolarPolarAngle = deg2rad(0); OptionsStruct.DipolarAzimuthAngle = 0; -OptionsStruct.ScatteringLength = 75.00; +OptionsStruct.ScatteringLength = 85; -OptionsStruct.TrapFrequencies = [0, 0, 500]; -OptionsStruct.TrapPotentialType = 'None'; +AspectRatio = 3.7; +HorizontalTrapFrequency = 125; +VerticalTrapFrequency = 125; +OptionsStruct.TrapFrequencies = [HorizontalTrapFrequency, HorizontalTrapFrequency, VerticalTrapFrequency]; +OptionsStruct.TrapPotentialType = 'Harmonic'; -OptionsStruct.NumberOfGridPoints = [128, 128]; -OptionsStruct.Dimensions = [lx, ly]; -OptionsStruct.TimeStepSize = 1E-3; % in s -OptionsStruct.MinimumTimeStepSize = 1E-5; % in s -OptionsStruct.TimeCutOff = 2E6; % in s +OptionsStruct.NumberOfGridPoints = [256, 256, 128]; +OptionsStruct.Dimensions = [30, 30, 15]; +OptionsStruct.CutoffType = 'Cylindrical'; +OptionsStruct.SimulationMode = 'ImaginaryTimeEvolution'; % 'ImaginaryTimeEvolution' | 'RealTimeEvolution' +OptionsStruct.TimeStepSize = 0.002; % in s +OptionsStruct.MinimumTimeStepSize = 1E-6; % in s +OptionsStruct.TimeCutOff = 1E6; % in s OptionsStruct.EnergyTolerance = 5E-10; OptionsStruct.ResidualTolerance = 1E-05; OptionsStruct.NoiseScaleFactor = 0.05; -OptionsStruct.IncludeDDICutOff = false; - -OptionsStruct.MaxIterations = 10; -OptionsStruct.VariationalWidth = 1.15; -OptionsStruct.WidthLowerBound = 0.01; -OptionsStruct.WidthUpperBound = 12; -OptionsStruct.WidthCutoff = 1e-2; OptionsStruct.PlotLive = false; -OptionsStruct.JobNumber = 1; +OptionsStruct.JobNumber = 0; OptionsStruct.RunOnGPU = true; OptionsStruct.SaveData = true; -OptionsStruct.SaveDirectory = './Results/Data_TiltingOfDipoles/AdjustedSystemSize/Hz500'; +OptionsStruct.SaveDirectory = sprintf('./Results/Data_3D/AspectRatio%s', strrep(num2str(AspectRatio), '.', '_')); options = Helper.convertstruct2cell(OptionsStruct); clear OptionsStruct -solver = VariationalSolver2D.DipolarGas(options{:}); -pot = VariationalSolver2D.Potentials(options{:}); -solver.Potential = pot.trap(); +sim = Simulator.DipolarGas(options{:}); +pot = Simulator.Potentials(options{:}); +sim.Potential = pot.trap(); % + pot.repulsive_chopstick(); -%-% Run Solver %-% -[Params, Transf, psi, V, VDk] = solver.run(); +%-% Run Simulation %-% +[Params, Transf, psi, V, VDk] = sim.run(); -%% v_z = 500, theta = 7.5: a_s = 75.00 - -a = 2.055; -scalingfactor = 5; -lx = scalingfactor*a; -ly = scalingfactor*sqrt(3)*a; - -% Initialize OptionsStruct -OptionsStruct = struct; - -% Assign values to OptionsStruct -OptionsStruct.NumberOfAtoms = ppum * (lx*ly); -OptionsStruct.DipolarPolarAngle = deg2rad(7.5); -OptionsStruct.DipolarAzimuthAngle = 0; -OptionsStruct.ScatteringLength = 75.00; - -OptionsStruct.TrapFrequencies = [0, 0, 500]; -OptionsStruct.TrapPotentialType = 'None'; - -OptionsStruct.NumberOfGridPoints = [128, 128]; -OptionsStruct.Dimensions = [lx, ly]; -OptionsStruct.TimeStepSize = 1E-3; % in s -OptionsStruct.MinimumTimeStepSize = 1E-5; % in s -OptionsStruct.TimeCutOff = 2E6; % in s -OptionsStruct.EnergyTolerance = 5E-10; -OptionsStruct.ResidualTolerance = 1E-05; -OptionsStruct.NoiseScaleFactor = 0.05; -OptionsStruct.IncludeDDICutOff = false; - -OptionsStruct.MaxIterations = 10; -OptionsStruct.VariationalWidth = 1.15; -OptionsStruct.WidthLowerBound = 0.01; -OptionsStruct.WidthUpperBound = 12; -OptionsStruct.WidthCutoff = 1e-2; - -OptionsStruct.PlotLive = false; -OptionsStruct.JobNumber = 2; -OptionsStruct.RunOnGPU = true; -OptionsStruct.SaveData = true; -OptionsStruct.SaveDirectory = './Results/Data_TiltingOfDipoles/AdjustedSystemSize/Hz500'; -options = Helper.convertstruct2cell(OptionsStruct); -clear OptionsStruct - -solver = VariationalSolver2D.DipolarGas(options{:}); -pot = VariationalSolver2D.Potentials(options{:}); -solver.Potential = pot.trap(); - -%-% Run Solver %-% -[Params, Transf, psi, V, VDk] = solver.run(); -%% v_z = 500, theta = 10: a_s = 75.00 - -a = 2.055; -scalingfactor = 5; -lx = scalingfactor*a; -ly = scalingfactor*sqrt(3)*a; - -% Initialize OptionsStruct -OptionsStruct = struct; - -% Assign values to OptionsStruct -OptionsStruct.NumberOfAtoms = ppum * (lx*ly); -OptionsStruct.DipolarPolarAngle = deg2rad(10); -OptionsStruct.DipolarAzimuthAngle = 0; -OptionsStruct.ScatteringLength = 75.00; - -OptionsStruct.TrapFrequencies = [0, 0, 500]; -OptionsStruct.TrapPotentialType = 'None'; - -OptionsStruct.NumberOfGridPoints = [128, 128]; -OptionsStruct.Dimensions = [lx, ly]; -OptionsStruct.TimeStepSize = 1E-3; % in s -OptionsStruct.MinimumTimeStepSize = 1E-5; % in s -OptionsStruct.TimeCutOff = 2E6; % in s -OptionsStruct.EnergyTolerance = 5E-10; -OptionsStruct.ResidualTolerance = 1E-05; -OptionsStruct.NoiseScaleFactor = 0.05; -OptionsStruct.IncludeDDICutOff = false; - -OptionsStruct.MaxIterations = 10; -OptionsStruct.VariationalWidth = 1.15; -OptionsStruct.WidthLowerBound = 0.01; -OptionsStruct.WidthUpperBound = 12; -OptionsStruct.WidthCutoff = 1e-2; - -OptionsStruct.PlotLive = false; -OptionsStruct.JobNumber = 3; -OptionsStruct.RunOnGPU = true; -OptionsStruct.SaveData = true; -OptionsStruct.SaveDirectory = './Results/Data_TiltingOfDipoles/AdjustedSystemSize/Hz500'; -options = Helper.convertstruct2cell(OptionsStruct); -clear OptionsStruct - -solver = VariationalSolver2D.DipolarGas(options{:}); -pot = VariationalSolver2D.Potentials(options{:}); -solver.Potential = pot.trap(); - -%-% Run Solver %-% -[Params, Transf, psi, V, VDk] = solver.run(); - -%% v_z = 500, theta = 15: a_s = 75.00 - -a = 2.0875; -scalingfactor = 5; -lx = scalingfactor*a; -ly = scalingfactor*sqrt(3)*a; - -% Initialize OptionsStruct -OptionsStruct = struct; - -% Assign values to OptionsStruct -OptionsStruct.NumberOfAtoms = ppum * (lx*ly); -OptionsStruct.DipolarPolarAngle = deg2rad(15); -OptionsStruct.DipolarAzimuthAngle = 0; -OptionsStruct.ScatteringLength = 75.00; - -OptionsStruct.TrapFrequencies = [0, 0, 500]; -OptionsStruct.TrapPotentialType = 'None'; - -OptionsStruct.NumberOfGridPoints = [128, 128]; -OptionsStruct.Dimensions = [lx, ly]; -OptionsStruct.TimeStepSize = 1E-3; % in s -OptionsStruct.MinimumTimeStepSize = 1E-5; % in s -OptionsStruct.TimeCutOff = 2E6; % in s -OptionsStruct.EnergyTolerance = 5E-10; -OptionsStruct.ResidualTolerance = 1E-05; -OptionsStruct.NoiseScaleFactor = 0.05; -OptionsStruct.IncludeDDICutOff = false; - -OptionsStruct.MaxIterations = 10; -OptionsStruct.VariationalWidth = 1.15; -OptionsStruct.WidthLowerBound = 0.01; -OptionsStruct.WidthUpperBound = 12; -OptionsStruct.WidthCutoff = 1e-2; - -OptionsStruct.PlotLive = false; -OptionsStruct.JobNumber = 4; -OptionsStruct.RunOnGPU = true; -OptionsStruct.SaveData = true; -OptionsStruct.SaveDirectory = './Results/Data_TiltingOfDipoles/AdjustedSystemSize/Hz500'; -options = Helper.convertstruct2cell(OptionsStruct); -clear OptionsStruct - -solver = VariationalSolver2D.DipolarGas(options{:}); -pot = VariationalSolver2D.Potentials(options{:}); -solver.Potential = pot.trap(); - -%-% Run Solver %-% -[Params, Transf, psi, V, VDk] = solver.run(); \ No newline at end of file diff --git a/Dipolar-Gas-Simulator/+Simulator/@Calculator/calculateChemicalPotential.m b/Dipolar-Gas-Simulator/+Simulator/@Calculator/calculateChemicalPotential.m index 96a6a84..6efd822 100644 --- a/Dipolar-Gas-Simulator/+Simulator/@Calculator/calculateChemicalPotential.m +++ b/Dipolar-Gas-Simulator/+Simulator/@Calculator/calculateChemicalPotential.m @@ -1,29 +1,29 @@ function muchem = calculateChemicalPotential(~,psi,Params,Transf,VDk,V) - %Parameters + % Parameters normfac = Params.Lx*Params.Ly*Params.Lz/numel(psi); - KEop= 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2); + KEop = 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2); % DDIs - frho=fftn(abs(psi).^2); - Phi=real(ifftn(frho.*VDk)); + frho = fftn(abs(psi).^2); + Phi = real(ifftn(frho.*VDk)); + Eddi = (Params.gdd*Phi.*abs(psi).^2); - Eddi = (Params.gdd*Phi.*abs(psi).^2); + % Kinetic energy + Ekin = KEop.*abs(fftn(psi)*normfac).^2; + Ekin = sum(Ekin(:))*Transf.dkx*Transf.dky*Transf.dkz/(2*pi)^3; - %Kinetic energy - Ekin = KEop.*abs(fftn(psi)*normfac).^2; - Ekin = trapz(Ekin(:))*Transf.dkx*Transf.dky*Transf.dkz/(2*pi)^3; + % Potential energy + Epot = V.*abs(psi).^2; - %Potential energy - Epot = V.*abs(psi).^2; + % Contact interactions + Eint = Params.gs*abs(psi).^4; - %Contact interactions - Eint = Params.gs*abs(psi).^4; + % Quantum fluctuations + Eqf = Params.gammaQF*abs(psi).^5; - %Quantum fluctuations - Eqf = Params.gammaQF*abs(psi).^5; - - %Total energy - muchem = Ekin + trapz(Epot(:) + Eint(:) + Eddi(:) + Eqf(:))*Transf.dx*Transf.dy*Transf.dz; % - muchem = muchem / Params.N; + % Total energy + muchem = Ekin + sum(Epot(:) + Eint(:) + Eddi(:) + Eqf(:))*Transf.dx*Transf.dy*Transf.dz; % + muchem = muchem / Params.N; + end \ No newline at end of file diff --git a/Dipolar-Gas-Simulator/+Simulator/@Calculator/calculateEnergyComponents.m b/Dipolar-Gas-Simulator/+Simulator/@Calculator/calculateEnergyComponents.m index 92dd32c..a8d53b6 100644 --- a/Dipolar-Gas-Simulator/+Simulator/@Calculator/calculateEnergyComponents.m +++ b/Dipolar-Gas-Simulator/+Simulator/@Calculator/calculateEnergyComponents.m @@ -1,35 +1,30 @@ function E = calculateEnergyComponents(~,psi,Params,Transf,VDk,V) - %Parameters - - KEop= 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2); + % Parameters + KEop = 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2); normfac = Params.Lx*Params.Ly*Params.Lz/numel(psi); % DDIs - frho = fftn(abs(psi).^2); - Phi = real(ifftn(frho.*VDk)); + frho = fftn(abs(psi).^2); + Phi = real(ifftn(frho.*VDk)); - Eddi = 0.5*Params.gdd*Phi.*abs(psi).^2; - E.Eddi = trapz(Eddi(:))*Transf.dx*Transf.dy*Transf.dz; + Eddi = 0.5*Params.gdd*Phi.*abs(psi).^2; + E.Eddi = sum(Eddi(:))*Transf.dx*Transf.dy*Transf.dz; - % EddiTot = trapz(Eddi(:))*Transf.dx*Transf.dy*Transf.dz; - - %Kinetic energy - % psik = ifftshift(fftn(fftshift(psi)))*normfac; - - Ekin = KEop.*abs(fftn(psi)*normfac).^2; - E.Ekin = trapz(Ekin(:))*Transf.dkx*Transf.dky*Transf.dkz/(2*pi)^3; + % Kinetic energy + Ekin = KEop.*abs(fftn(psi)*normfac).^2; + E.Ekin = sum(Ekin(:))*Transf.dkx*Transf.dky*Transf.dkz/(2*pi)^3; % Potential energy - Epot = V.*abs(psi).^2; - E.Epot = trapz(Epot(:))*Transf.dx*Transf.dy*Transf.dz; + Epot = V.*abs(psi).^2; + E.Epot = sum(Epot(:))*Transf.dx*Transf.dy*Transf.dz; %Contact interactions - Eint = 0.5*Params.gs*abs(psi).^4; - E.Eint = trapz(Eint(:))*Transf.dx*Transf.dy*Transf.dz; + Eint = 0.5*Params.gs*abs(psi).^4; + E.Eint = sum(Eint(:))*Transf.dx*Transf.dy*Transf.dz; %Quantum fluctuations - Eqf = 0.4*Params.gammaQF*abs(psi).^5; - E.Eqf = trapz(Eqf(:))*Transf.dx*Transf.dy*Transf.dz; + Eqf = 0.4*Params.gammaQF*abs(psi).^5; + E.Eqf = sum(Eqf(:))*Transf.dx*Transf.dy*Transf.dz; end diff --git a/Dipolar-Gas-Simulator/+Simulator/@Calculator/calculateNormalizedResiduals.m b/Dipolar-Gas-Simulator/+Simulator/@Calculator/calculateNormalizedResiduals.m index 134e426..b1f2fa3 100644 --- a/Dipolar-Gas-Simulator/+Simulator/@Calculator/calculateNormalizedResiduals.m +++ b/Dipolar-Gas-Simulator/+Simulator/@Calculator/calculateNormalizedResiduals.m @@ -1,25 +1,26 @@ function res = calculateNormalizedResiduals(~,psi,Params,Transf,VDk,V,muchem) - KEop= 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2); + KEop = 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2); % DDIs - frho=fftn(abs(psi).^2); - Phi=real(ifftn(frho.*VDk)); + frho = fftn(abs(psi).^2); + Phi = real(ifftn(frho.*VDk)); Eddi = Params.gdd*Phi.*psi; - %Kinetic energy + % Kinetic energy Ekin = ifftn(KEop.*fftn(psi)); - %Potential energy + % Potential energy Epot = V.*psi; - %Contact interactions + % Contact interactions Eint = Params.gs*abs(psi).^2.*psi; - %Quantum fluctuations - Eqf = Params.gammaQF*abs(psi).^3.*psi; + % Quantum fluctuations + Eqf = Params.gammaQF*abs(psi).^3.*psi; - %Total energy - res = trapz(abs(Ekin(:) + Epot(:) + Eint(:) + Eddi(:) + Eqf(:) - muchem*psi(:))*Transf.dx*Transf.dy*Transf.dz)/trapz(abs(muchem*psi(:))*Transf.dx*Transf.dy*Transf.dz); + % Total energy + res = sum(abs(Ekin(:) + Epot(:) + Eint(:) + Eddi(:) + Eqf(:) - muchem*psi(:))*Transf.dx*Transf.dy*Transf.dz)/sum(abs(muchem*psi(:))*Transf.dx*Transf.dy*Transf.dz); + end \ No newline at end of file diff --git a/Dipolar-Gas-Simulator/+Simulator/@Calculator/calculateNumericalHankelTransform.m b/Dipolar-Gas-Simulator/+Simulator/@Calculator/calculateNumericalHankelTransform.m index b3d96c6..d3d6550 100644 --- a/Dipolar-Gas-Simulator/+Simulator/@Calculator/calculateNumericalHankelTransform.m +++ b/Dipolar-Gas-Simulator/+Simulator/@Calculator/calculateNumericalHankelTransform.m @@ -5,25 +5,26 @@ function VDkSemi = calculateNumericalHankelTransform(~,kr,kz,Rmax,Zmax,Nr) Nr = 5e4; end - Nz = 64; + Nz = 64; farRmultiple = 2000; % midpoint grids for the integration over 0 0))); - addParameter(p, 'NumberOfTimeSteps', 2e6,... + addParameter(p, 'MinimumTimeStepSize', 1e-6,... + @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); + addParameter(p, 'TimeCutOff', 2e6,... @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); addParameter(p, 'EnergyTolerance', 1e-10,... @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); addParameter(p, 'ResidualTolerance', 1e-10,... - @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); - addParameter(p, 'MinimumTimeStepSize', 1e-6,... @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); + addParameter(p, 'NoiseScaleFactor', 4,... + @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); + + addParameter(p, 'IncludeDDICutOff', true,... + @islogical); + addParameter(p, 'PlotLive', false,... + @islogical); addParameter(p, 'JobNumber', 1,... @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); addParameter(p, 'RunOnGPU', false,... @@ -87,11 +96,14 @@ classdef DipolarGas < handle & matlab.mixin.Copyable this.Potential = NaN; this.SimulationMode = p.Results.SimulationMode; this.TimeStepSize = p.Results.TimeStepSize; - this.NumberOfTimeSteps = p.Results.NumberOfTimeSteps; + this.MinimumTimeStepSize = p.Results.MinimumTimeStepSize; + this.TimeCutOff = p.Results.TimeCutOff; this.EnergyTolerance = p.Results.EnergyTolerance; this.ResidualTolerance = p.Results.ResidualTolerance; - this.MinimumTimeStepSize = p.Results.MinimumTimeStepSize; + this.NoiseScaleFactor = p.Results.NoiseScaleFactor; + this.IncludeDDICutOff = p.Results.IncludeDDICutOff; + this.PlotLive = p.Results.PlotLive; this.JobNumber = p.Results.JobNumber; this.RunOnGPU = p.Results.RunOnGPU; this.DebugMode = p.Results.DebugMode; @@ -113,12 +125,12 @@ classdef DipolarGas < handle & matlab.mixin.Copyable function ret = get.TimeStepSize(this) ret = this.TimeStepSize; end - function set.NumberOfTimeSteps(this, val) + function set.TimeCutOff(this, val) assert(val <= 2E6, 'Not time efficient to compute for time spans longer than 2E6 seconds!'); - this.NumberOfTimeSteps = val; + this.TimeCutOff = val; end - function ret = get.NumberOfTimeSteps(this) - ret = this.NumberOfTimeSteps; + function ret = get.TimeCutOff(this) + ret = this.TimeCutOff; end function set.NumberOfAtoms(this, val) assert(val <= 1E6, '!!Not time efficient to compute for atom numbers larger than 1,000,000!!'); diff --git a/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/initialize.m b/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/initialize.m index 2731a5b..01f3cd9 100644 --- a/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/initialize.m +++ b/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/initialize.m @@ -8,8 +8,12 @@ function [psi,V,VDk] = initialize(this,Params,Transf,TransfRad) if isfile(strcat(this.SaveDirectory, '/VDk_M.mat')) VDk = load(sprintf(strcat(this.SaveDirectory, '/VDk_M.mat'))); VDk = VDk.VDk; + if ~isequal(size(VDk), this.NumberOfGridPoints) + VDk = this.Calculator.calculateVDk(Params,Transf,TransfRad, this.IncludeDDICutOff); + save(sprintf(strcat(this.SaveDirectory, '/VDk_M.mat')),'VDk'); + end else - VDk = this.Calculator.calculateVDkWithCutoff(Params,Transf,TransfRad); + VDk = this.Calculator.calculateVDk(Params,Transf,TransfRad, this.IncludeDDICutOff); save(sprintf(strcat(this.SaveDirectory, '/VDk_M.mat')),'VDk'); end fprintf('Computed and saved DDI potential in Fourier space with %s cutoff.\n', this.Calculator.CutoffType) diff --git a/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/propagateWavefunction.m b/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/propagateWavefunction.m index 7bad664..d1ee687 100644 --- a/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/propagateWavefunction.m +++ b/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/propagateWavefunction.m @@ -4,63 +4,70 @@ function [psi] = propagateWavefunction(this,psi,Params,Transf,VDk,V,t_idx,Observ switch this.SimulationMode case 'ImaginaryTimeEvolution' - dt=-1j*abs(this.TimeStepSize); + dt = -1j*abs(this.TimeStepSize); % Imaginary Time - KEop= 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2); - Observ.residual = 1; Observ.res = 1; + KEop = 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2); + Observ.residual = 1; + Observ.res = 1; - muchem = this.Calculator.calculateChemicalPotential(psi,Params,Transf,VDk,V); - AdaptIdx = 0; + muchem = this.Calculator.calculateChemicalPotential(psi,Params,Transf,VDk,V); + AdaptIdx = 0; - pb = Helper.ProgressBar(); - pb.run('Running simulation in imaginary time: '); + if this.PlotLive + Plotter.plotLive(psi,Params,Transf,Observ) + drawnow + end while t_idx < Params.sim_time_cut_off - %kin - psi = fftn(psi); - psi = psi.*exp(-0.5*1i*dt*KEop); - psi = ifftn(psi); + % kin + psi = fftn(psi); + psi = psi.*exp(-0.5*1i*dt*KEop); + psi = ifftn(psi); - %DDI - frho = fftn(abs(psi).^2); - Phi = real(ifftn(frho.*VDk)); + % DDI + frho = fftn(abs(psi).^2); + Phi = real(ifftn(frho.*VDk)); - %Real-space - psi = psi.*exp(-1i*dt*(V + Params.gs*abs(psi).^2 + Params.gdd*Phi + Params.gammaQF*abs(psi).^3 - muchem)); + % Real-space + psi = psi.*exp(-1i*dt*(V + Params.gs*abs(psi).^2 + Params.gdd*Phi + Params.gammaQF*abs(psi).^3 - muchem)); - %kin - psi = fftn(psi); - psi = psi.*exp(-0.5*1i*dt*KEop); - psi = ifftn(psi); + % kin + psi = fftn(psi); + psi = psi.*exp(-0.5*1i*dt*KEop); + psi = ifftn(psi); - %Renorm - Norm = trapz(abs(psi(:)).^2)*Transf.dx*Transf.dy*Transf.dz; - psi = sqrt(Params.N)*psi/sqrt(Norm); + % Renorm + Norm = sum(abs(psi(:)).^2)*Transf.dx*Transf.dy*Transf.dz; + psi = sqrt(Params.N)*psi/sqrt(Norm); - muchem = this.Calculator.calculateChemicalPotential(psi,Params,Transf,VDk,V); + muchem = this.Calculator.calculateChemicalPotential(psi,Params,Transf,VDk,V); - if mod(t_idx,1000) == 0 + if mod(t_idx,500) == 0 - %Change in Energy - E = this.Calculator.calculateTotalEnergy(psi,Params,Transf,VDk,V); - E = E/Norm; - Observ.EVec = [Observ.EVec E]; + % 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]; + % Chemical potential + Observ.mucVec = [Observ.mucVec muchem]; - %Normalized residuals - res = this.Calculator.calculateNormalizedResiduals(psi,Params,Transf,VDk,V,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; - Observ.res_idx = Observ.res_idx + 1; - - save(sprintf('./Data/Run_%03i/psi_gs.mat',Params.njob),'psi','muchem','Observ','t_idx','Transf','Params','VDk','V'); + 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','t_idx','Transf','Params','VDk','V'); %Adaptive time step -- Careful, this can quickly get out of control relres = abs(Observ.residual(Observ.res_idx)-Observ.residual(Observ.res_idx-1))/Observ.residual(Observ.res_idx); - if relres <1e-5 + if relres < 1e-5 if AdaptIdx > 4 && abs(dt) > Params.mindt dt = dt / 2; fprintf('Time step changed to '); disp(dt); @@ -79,23 +86,22 @@ function [psi] = propagateWavefunction(this,psi,Params,Transf,VDk,V,t_idx,Observ break end t_idx=t_idx+1; - pb.run(100*t_idx/Params.sim_time_cut_off); end - %Change in Energy - E = this.Calculator.calculateTotalEnergy(psi,Params,Transf,VDk,V); - E = E/Norm; - Observ.EVec = [Observ.EVec E]; + % Change in Energy + E = this.Calculator.calculateTotalEnergy(psi,Params,Transf,VDk,V); + E = E/Norm; + Observ.EVec = [Observ.EVec E]; % Phase coherence - [PhaseC] = this.Calculator.calculatePhaseCoherence(psi,Transf,Params); - Observ.PCVec = [Observ.PCVec PhaseC]; + [PhaseC] = this.Calculator.calculatePhaseCoherence(psi,Transf,Params); + Observ.PCVec = [Observ.PCVec PhaseC]; Observ.res_idx = Observ.res_idx + 1; - pb.run(' - Job Completed!'); + disp('Run completed!'); disp('Saving data...'); - save(sprintf('./Data/Run_%03i/psi_gs.mat',Params.njob),'psi','muchem','Observ','t_idx','Transf','Params','VDk','V'); + save(sprintf(strcat(this.SaveDirectory, '/Run_%03i/psi_gs.mat'),Params.njob),'psi','muchem','Observ','t_idx','Transf','Params','VDk','V'); disp('Save complete!'); case 'RealTimeEvolution' @@ -106,9 +112,6 @@ function [psi] = propagateWavefunction(this,psi,Params,Transf,VDk,V,t_idx,Observ muchem = this.Calculator.calculateChemicalPotential(psi,Params,Transf,VDk,V); - pb = Helper.ProgressBar(); - pb.run('Running simulation in real time: '); - while t_idx < Params.sim_time_cut_off % Parameters at time t @@ -134,7 +137,7 @@ function [psi] = propagateWavefunction(this,psi,Params,Transf,VDk,V,t_idx,Observ if mod(t_idx,1000)==0 %Change in Normalization - Norm = trapz(abs(psi(:)).^2)*Transf.dx*Transf.dy*Transf.dz; % normalisation + Norm = sum(abs(psi(:)).^2)*Transf.dx*Transf.dy*Transf.dz; % normalisation Observ.NormVec = [Observ.NormVec Norm]; %Change in Energy @@ -156,11 +159,10 @@ function [psi] = propagateWavefunction(this,psi,Params,Transf,VDk,V,t_idx,Observ break end t_idx=t_idx+1; - pb.run(100*t_idx/Params.sim_time_cut_off); end %Change in Normalization - Norm = trapz(abs(psi(:)).^2)*Transf.dx*Transf.dy*Transf.dz; % normalisation + Norm = sum(abs(psi(:)).^2)*Transf.dx*Transf.dy*Transf.dz; % normalisation Observ.NormVec = [Observ.NormVec Norm]; %Change in Energy @@ -175,7 +177,7 @@ function [psi] = propagateWavefunction(this,psi,Params,Transf,VDk,V,t_idx,Observ Observ.tVecPlot = [Observ.tVecPlot tVal]; Observ.res_idx = Observ.res_idx + 1; - pb.run(' - Run Completed!'); + disp('Run completed!'); disp('Saving data...'); save(sprintf('./Data/Run_%03i/TimeEvolution/psi_%i.mat',Params.njob,Observ.res_idx),'psi','muchem','Observ','t_idx'); disp('Save complete!'); diff --git a/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/run.m b/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/run.m index 1d4adfb..3f4dcf9 100644 --- a/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/run.m +++ b/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/run.m @@ -1,22 +1,23 @@ function [Params, Transf, psi,V,VDk] = run(this) % --- Obtain simulation parameters --- - [Params] = this.setupParameters(); + [Params] = this.setupParameters(); this.SimulationParameters = Params; % --- Set up spatial grids and transforms --- - [Transf] = this.setupSpace(Params); - [TransfRad] = this.setupSpaceRadial(Params); + [Transf] = this.setupSpace(Params); + [TransfRad] = this.setupSpaceRadial(Params); % --- Initialize --- mkdir(sprintf(this.SaveDirectory)) - [psi,V,VDk] = this.initialize(Params,Transf,TransfRad); + [psi,V,VDk] = this.initialize(Params,Transf,TransfRad); Observ.EVec = []; Observ.NormVec = []; Observ.PCVec = []; Observ.tVecPlot = []; Observ.mucVec = []; - t_idx = 1; % Start at t = 0; - Observ.res_idx = 1; + t_idx = 1; % Start at t = 0; + Observ.res_idx = 1; % --- Run Simulation --- - mkdir(sprintf('./Data/Run_%03i',Params.njob)) + mkdir(sprintf(this.SaveDirectory)) + mkdir(sprintf(strcat(this.SaveDirectory, '/Run_%03i'),Params.njob)) [psi] = this.propagateWavefunction(psi,Params,Transf,VDk,V,t_idx,Observ); end diff --git a/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/setupParameters.m b/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/setupParameters.m index cf41fe0..0487153 100644 --- a/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/setupParameters.m +++ b/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/setupParameters.m @@ -6,7 +6,7 @@ function [Params] = setupParameters(this) muB = CONSTANTS.BohrMagneton; % [J/T] a0 = CONSTANTS.BohrRadius; % [m] m0 = CONSTANTS.AtomicMassUnit; % [kg] - w0 = 2*pi*100; % Angular frequency unit [s^-1] + w0 = 2*pi*61.658214297935530; % Angular frequency unit [s^-1] mu0factor = 0.3049584233607396; % =(m0/me)*pi*alpha^2 -- me=mass of electron, alpha=fine struct. const. % mu0=mu0factor *hbar^2*a0/(m0*muB^2) % Number of points in each direction @@ -21,7 +21,6 @@ function [Params] = setupParameters(this) % Mass, length scale Params.m = CONSTANTS.Dy164Mass; - l0 = sqrt(hbar/(Params.m*w0)); % Defining a harmonic oscillator length % Atom numbers Params.N = this.NumberOfAtoms; @@ -44,14 +43,18 @@ function [Params] = setupParameters(this) % Tolerances Params.Etol = this.EnergyTolerance; Params.rtol = this.ResidualTolerance; - Params.sim_time_cut_off = this.NumberOfTimeSteps; % sometimes the imaginary time gets a little stuck - % even though the solution is good, this just stops it going on forever - Params.mindt = this.MinimumTimeStepSize; % Minimum size for a time step using adaptive dt + Params.sim_time_cut_off = this.TimeCutOff; % sometimes the imaginary time gets a little stuck + % even though the solution is good, this just stops it going on forever + Params.mindt = this.MinimumTimeStepSize; % Minimum size for a time step using adaptive dt - Params.njob = this.JobNumber; + Params.nsf = this.NoiseScaleFactor; + + Params.njob = this.JobNumber; % ================ Parameters defined by those above ================ % - + % Length scale + l0 = sqrt(hbar/(Params.m*w0)); % Defining a harmonic oscillator length + % Contact interaction strength (units of l0/m) Params.gs = 4*pi*Params.as/l0; @@ -59,7 +62,7 @@ function [Params] = setupParameters(this) Params.add = mu0*Params.mu^2*Params.m/(12*pi*hbar^2); % DDI strength - Params.gdd = 12*pi*Params.add/l0; %sometimes the 12 is a 4 --> depends on how Vdk (DDI) is defined + Params.gdd = 4*pi*Params.add/l0; %sometimes the 12 is a 4 --> depends on how Vdk (DDI) is defined % Trap gamma Params.gx = (Params.wx/w0)^2; diff --git a/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/setupSpace.m b/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/setupSpace.m index bf58f64..9ce6d8f 100644 --- a/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/setupSpace.m +++ b/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/setupSpace.m @@ -1,33 +1,48 @@ function [Transf] = setupSpace(~,Params) - Transf.Xmax = 0.5*Params.Lx; - Transf.Ymax = 0.5*Params.Ly; - Transf.Zmax = 0.5*Params.Lz; - Nz = Params.Nz; Nx = Params.Nx; Ny = Params.Ny; + Transf.Xmax = 0.5*Params.Lx; + Transf.Ymax = 0.5*Params.Ly; + Transf.Zmax = 0.5*Params.Lz; + + Nz = Params.Nz; + Nx = Params.Nx; + Ny = Params.Ny; % Fourier grids - x = linspace(-0.5*Params.Lx,0.5*Params.Lx-Params.Lx/Params.Nx,Params.Nx); - Kmax = pi*Params.Nx/Params.Lx; - kx = linspace(-Kmax,Kmax,Nx+1); - kx = kx(1:end-1); dkx = kx(2)-kx(1); - kx = fftshift(kx); + x = linspace(-0.5*Params.Lx,0.5*Params.Lx-Params.Lx/Params.Nx,Params.Nx); + Kmax = pi*Params.Nx/Params.Lx; + kx = linspace(-Kmax,Kmax,Nx+1); + kx = kx(1:end-1); + dkx = kx(2)-kx(1); + kx = fftshift(kx); - y = linspace(-0.5*Params.Ly,0.5*Params.Ly-Params.Ly/Params.Ny,Params.Ny); - Kmax = pi*Params.Ny/Params.Ly; - ky = linspace(-Kmax,Kmax,Ny+1); - ky = ky(1:end-1); dky = ky(2)-ky(1); - ky = fftshift(ky); + y = linspace(-0.5*Params.Ly,0.5*Params.Ly-Params.Ly/Params.Ny,Params.Ny); + Kmax = pi*Params.Ny/Params.Ly; + ky = linspace(-Kmax,Kmax,Ny+1); + ky = ky(1:end-1); + dky = ky(2)-ky(1); + ky = fftshift(ky); - z = linspace(-0.5*Params.Lz,0.5*Params.Lz-Params.Lz/Params.Nz,Params.Nz); - Kmax = pi*Params.Nz/Params.Lz; - kz = linspace(-Kmax,Kmax,Nz+1); - kz = kz(1:end-1); dkz = kz(2)-kz(1); - kz = fftshift(kz); + z = linspace(-0.5*Params.Lz,0.5*Params.Lz-Params.Lz/Params.Nz,Params.Nz); + Kmax = pi*Params.Nz/Params.Lz; + kz = linspace(-Kmax,Kmax,Nz+1); + kz = kz(1:end-1); + dkz = kz(2)-kz(1); + kz = fftshift(kz); + + [Transf.X,Transf.Y,Transf.Z] = ndgrid(x,y,z); + [Transf.KX,Transf.KY,Transf.KZ] = ndgrid(kx,ky,kz); + Transf.x = x; + Transf.y = y; + Transf.z = z; + Transf.kx = kx; + Transf.ky = ky; + Transf.kz = kz; + Transf.dx = x(2)-x(1); + Transf.dy = y(2)-y(1); + Transf.dz = z(2)-z(1); + Transf.dkx = dkx; + Transf.dky = dky; + Transf.dkz = dkz; - [Transf.X,Transf.Y,Transf.Z]=ndgrid(x,y,z); - [Transf.KX,Transf.KY,Transf.KZ]=ndgrid(kx,ky,kz); - Transf.x = x; Transf.y = y; Transf.z = z; - Transf.kx = kx; Transf.ky = ky; Transf.kz = kz; - Transf.dx = x(2)-x(1); Transf.dy = y(2)-y(1); Transf.dz = z(2)-z(1); - Transf.dkx = dkx; Transf.dky = dky; Transf.dkz = dkz; 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 3d604a4..15ec981 100644 --- a/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/setupWavefunction.m +++ b/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/setupWavefunction.m @@ -1,25 +1,32 @@ function [psi] = setupWavefunction(~,Params,Transf) - X = Transf.X; Y = Transf.Y; Z = Transf.Z; + 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; - Rx = 4*ellx; Ry = 4*elly; Rz = 4*ellz; - X0 = 0.0*Transf.Xmax; Y0 = 0.0*Transf.Ymax; Z0 = 0*Transf.Zmax; + Rx = 8*ellx; + Ry = 8*elly; + Rz = 8*ellz; + X0 = 0.0*Transf.Xmax; + Y0 = 0.0*Transf.Ymax; + Z0 = 0*Transf.Zmax; psi = exp(-(X-X0).^2/Rx^2-(Y-Y0).^2/Ry^2-(Z-Z0).^2/Rz^2); - cur_norm = trapz(abs(psi(:)).^2)*Transf.dx*Transf.dy*Transf.dz; + cur_norm = sum(abs(psi(:)).^2)*Transf.dx*Transf.dy*Transf.dz; psi = psi/sqrt(cur_norm); - % add some noise + % --- Adding some noise --- r = normrnd(0,1,size(X)); theta = rand(size(X)); noise = r.*exp(2*pi*1i*theta); - psi = psi + 0.01*noise; + psi = psi + Params.nsf*noise; - Norm = trapz(abs(psi(:)).^2)*Transf.dx*Transf.dy*Transf.dz; - psi = sqrt(Params.N)*psi/sqrt(Norm); + % Renormalize wavefunction + 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