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

Added custom cylindrical cut off, yet to be fully tested, minor mods to plotting routines.

Browse Source
This commit is contained in:
Karthik 2025-03-26 19:50:59 +01:00
parent db804b0117
commit 88f2554b8e
7 changed files with 51 additions and 24 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
Dipolar-Gas-Simulator/+Plotter/plotLive.m
View File

@ -32,7 +32,7 @@ function plotLive(psi,Params,Transf,Observ)
colormap(gca, Helper.Colormaps.plasma()) colormap(gca, Helper.Colormaps.plasma())
xlabel('$x$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 14) xlabel('$x$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 14)
ylabel('$z$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 14) ylabel('$z$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 14)
title('$|\Psi(x,y)|^2$', 'Interpreter', 'latex', 'FontSize', 14) title('$|\Psi(x,z)|^2$', 'Interpreter', 'latex', 'FontSize', 14)
nexttile; nexttile;
plotyz = pcolor(y,z,nyz'); plotyz = pcolor(y,z,nyz');
@ -43,7 +43,7 @@ function plotLive(psi,Params,Transf,Observ)
colormap(gca, Helper.Colormaps.plasma()) colormap(gca, Helper.Colormaps.plasma())
xlabel('$y$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 14) xlabel('$y$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 14)
ylabel('$z$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 14) ylabel('$z$ ($\mu$m)', 'Interpreter', 'latex', 'FontSize', 14)
title('$|\Psi(x,y)|^2$', 'Interpreter', 'latex', 'FontSize', 14) title('$|\Psi(y,z)|^2$', 'Interpreter', 'latex', 'FontSize', 14)
nexttile; nexttile;
plotxy = pcolor(x,y,nxy'); plotxy = pcolor(x,y,nxy');

6
Dipolar-Gas-Simulator/+Plotter/visualizeGSWavefunction.m
View File

@ -1,9 +1,9 @@
function visualizeGSWavefunction(run_index) function visualizeGSWavefunction(folder_path, run_index)
set(0,'defaulttextInterpreter','latex') set(0,'defaulttextInterpreter','latex')
set(groot, 'defaultAxesTickLabelInterpreter','latex'); set(groot, 'defaultLegendInterpreter','latex'); set(groot, 'defaultAxesTickLabelInterpreter','latex'); set(groot, 'defaultLegendInterpreter','latex');
folder_path = './Data'; % Load data
Data = load(sprintf(horzcat(folder_path, '/Run_%03i/psi_gs.mat'),run_index),'psi','Params','Transf','Observ'); Data = load(sprintf(horzcat(folder_path, '/Run_%03i/psi_gs.mat'),run_index),'psi','Params','Transf','Observ');
Params = Data.Params; Params = Data.Params;

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

@ -11,14 +11,17 @@ OptionsStruct.DipolarPolarAngle = deg2rad(0);
OptionsStruct.DipolarAzimuthAngle = 0; OptionsStruct.DipolarAzimuthAngle = 0;
OptionsStruct.ScatteringLength = 85; OptionsStruct.ScatteringLength = 85;
OptionsStruct.TrapFrequencies = [125, 125, 350]; AspectRatio = 2.8;
HorizontalTrapFrequency = 125;
VerticalTrapFrequency = AspectRatio * HorizontalTrapFrequency;
OptionsStruct.TrapFrequencies = [HorizontalTrapFrequency, HorizontalTrapFrequency, VerticalTrapFrequency];
OptionsStruct.TrapPotentialType = 'Harmonic'; OptionsStruct.TrapPotentialType = 'Harmonic';
OptionsStruct.NumberOfGridPoints = [128, 128, 64]; OptionsStruct.NumberOfGridPoints = [256, 256, 128];
OptionsStruct.Dimensions = [30, 30, 15]; OptionsStruct.Dimensions = [30, 30, 15];
OptionsStruct.CutoffType = 'Cylindrical'; OptionsStruct.CutoffType = 'Cylindrical';
OptionsStruct.SimulationMode = 'ImaginaryTimeEvolution'; % 'ImaginaryTimeEvolution' | 'RealTimeEvolution' OptionsStruct.SimulationMode = 'ImaginaryTimeEvolution'; % 'ImaginaryTimeEvolution' | 'RealTimeEvolution'
OptionsStruct.TimeStepSize = 0.002; % in s OptionsStruct.TimeStepSize = 0.001; % in s
OptionsStruct.MinimumTimeStepSize = 1E-6; % in s OptionsStruct.MinimumTimeStepSize = 1E-6; % in s
OptionsStruct.TimeCutOff = 1E6; % in s OptionsStruct.TimeCutOff = 1E6; % in s
OptionsStruct.EnergyTolerance = 5E-10; OptionsStruct.EnergyTolerance = 5E-10;
@ -26,10 +29,10 @@ OptionsStruct.ResidualTolerance = 1E-05;
OptionsStruct.NoiseScaleFactor = 0.05; OptionsStruct.NoiseScaleFactor = 0.05;
OptionsStruct.PlotLive = true; OptionsStruct.PlotLive = true;
OptionsStruct.JobNumber = 1; OptionsStruct.JobNumber = 0;
OptionsStruct.RunOnGPU = false; OptionsStruct.RunOnGPU = false;
OptionsStruct.SaveData = true; OptionsStruct.SaveData = true;
OptionsStruct.SaveDirectory = './Results/Data_3D'; OptionsStruct.SaveDirectory = sprintf('./Results/Data_3D/AspectRatio%s', strrep(num2str(AspectRatio), '.', '_'));
options = Helper.convertstruct2cell(OptionsStruct); options = Helper.convertstruct2cell(OptionsStruct);
clear OptionsStruct clear OptionsStruct
@ -47,7 +50,9 @@ 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 %% - Plot GS wavefunction
Plotter.visualizeGSWavefunction(Params.njob) SaveDirectory = './Results/Data_3D/AspectRatio2_8';
JobNumber = 0;
Plotter.visualizeGSWavefunction(SaveDirectory, JobNumber)
%% %%

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

@ -17,7 +17,7 @@ OptionsStruct.NumberOfGridPoints = [256, 256, 128];
OptionsStruct.Dimensions = [30, 30, 15]; OptionsStruct.Dimensions = [30, 30, 15];
OptionsStruct.CutoffType = 'Cylindrical'; OptionsStruct.CutoffType = 'Cylindrical';
OptionsStruct.SimulationMode = 'ImaginaryTimeEvolution'; % 'ImaginaryTimeEvolution' | 'RealTimeEvolution' OptionsStruct.SimulationMode = 'ImaginaryTimeEvolution'; % 'ImaginaryTimeEvolution' | 'RealTimeEvolution'
OptionsStruct.TimeStepSize = 0.002; % in s OptionsStruct.TimeStepSize = 0.001; % in s
OptionsStruct.MinimumTimeStepSize = 1E-6; % in s OptionsStruct.MinimumTimeStepSize = 1E-6; % in s
OptionsStruct.TimeCutOff = 1E6; % in s OptionsStruct.TimeCutOff = 1E6; % in s
OptionsStruct.EnergyTolerance = 5E-10; OptionsStruct.EnergyTolerance = 5E-10;
@ -58,7 +58,7 @@ OptionsStruct.NumberOfGridPoints = [256, 256, 128];
OptionsStruct.Dimensions = [30, 30, 15]; OptionsStruct.Dimensions = [30, 30, 15];
OptionsStruct.CutoffType = 'Cylindrical'; OptionsStruct.CutoffType = 'Cylindrical';
OptionsStruct.SimulationMode = 'ImaginaryTimeEvolution'; % 'ImaginaryTimeEvolution' | 'RealTimeEvolution' OptionsStruct.SimulationMode = 'ImaginaryTimeEvolution'; % 'ImaginaryTimeEvolution' | 'RealTimeEvolution'
OptionsStruct.TimeStepSize = 0.002; % in s OptionsStruct.TimeStepSize = 0.001; % in s
OptionsStruct.MinimumTimeStepSize = 1E-6; % in s OptionsStruct.MinimumTimeStepSize = 1E-6; % in s
OptionsStruct.TimeCutOff = 1E6; % in s OptionsStruct.TimeCutOff = 1E6; % in s
OptionsStruct.EnergyTolerance = 5E-10; OptionsStruct.EnergyTolerance = 5E-10;
@ -66,7 +66,7 @@ OptionsStruct.ResidualTolerance = 1E-05;
OptionsStruct.NoiseScaleFactor = 0.05; OptionsStruct.NoiseScaleFactor = 0.05;
OptionsStruct.PlotLive = false; OptionsStruct.PlotLive = false;
OptionsStruct.JobNumber = 1; OptionsStruct.JobNumber = 0;
OptionsStruct.RunOnGPU = true; OptionsStruct.RunOnGPU = true;
OptionsStruct.SaveData = true; OptionsStruct.SaveData = true;
OptionsStruct.SaveDirectory = sprintf('./Results/Data_3D/AspectRatio%s', strrep(num2str(AspectRatio), '.', '_')); OptionsStruct.SaveDirectory = sprintf('./Results/Data_3D/AspectRatio%s', strrep(num2str(AspectRatio), '.', '_'));

26
Dipolar-Gas-Simulator/+Simulator/@Calculator/Calculator.m
View File

@ -8,7 +8,9 @@ classdef Calculator < handle & matlab.mixin.Copyable
'OrderParameter', 1, ... 'OrderParameter', 1, ...
'PhaseCoherence', 1, ... 'PhaseCoherence', 1, ...
'TotalEnergy', 1, ... 'TotalEnergy', 1, ...
'CutoffType', 'Cylindrical'); 'CutoffType', 'Cylindrical', ...
'CustomCylindricalCutOffRadius', 1, ...
'CustomCylindricalCutOffHeight', 1);
end end
@ -20,6 +22,8 @@ classdef Calculator < handle & matlab.mixin.Copyable
PhaseCoherence; PhaseCoherence;
TotalEnergy; TotalEnergy;
CutoffType; CutoffType;
CustomCylindricalCutOffRadius
CustomCylindricalCutOffHeight
end end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@ -31,7 +35,7 @@ classdef Calculator < handle & matlab.mixin.Copyable
p = inputParser; p = inputParser;
p.KeepUnmatched = true; p.KeepUnmatched = true;
addParameter(p, 'CutoffType', this.CalculatorDefaults.CutoffType,... addParameter(p, 'CutoffType', this.CalculatorDefaults.CutoffType,...
@(x) any(strcmpi(x,{'Cylindrical','CylindricalInfiniteZ', 'Spherical'}))); @(x) any(strcmpi(x,{'Cylindrical','CylindricalInfiniteZ', 'Spherical', 'CustomCylindrical'})));
p.parse(varargin{:}); p.parse(varargin{:});
@ -42,16 +46,20 @@ classdef Calculator < handle & matlab.mixin.Copyable
this.PhaseCoherence = this.CalculatorDefaults.PhaseCoherence; this.PhaseCoherence = this.CalculatorDefaults.PhaseCoherence;
this.TotalEnergy = this.CalculatorDefaults.TotalEnergy; this.TotalEnergy = this.CalculatorDefaults.TotalEnergy;
this.CutoffType = p.Results.CutoffType; this.CutoffType = p.Results.CutoffType;
this.CustomCylindricalCutOffRadius = this.CalculatorDefaults.CustomCylindricalCutOffRadius;
this.CustomCylindricalCutOffHeight = this.CalculatorDefaults.CustomCylindricalCutOffHeight;
end end
function restoreDefaults(this) function restoreDefaults(this)
this.ChemicalPotential = this.CalculatorDefaults.ChemicalPotential; this.ChemicalPotential = this.CalculatorDefaults.ChemicalPotential;
this.EnergyComponents = this.CalculatorDefaults.EnergyComponents; this.EnergyComponents = this.CalculatorDefaults.EnergyComponents;
this.NormalizedResiduals = this.CalculatorDefaults.NormalizedResiduals; this.NormalizedResiduals = this.CalculatorDefaults.NormalizedResiduals;
this.OrderParameter = this.CalculatorDefaults.OrderParameter; this.OrderParameter = this.CalculatorDefaults.OrderParameter;
this.PhaseCoherence = this.CalculatorDefaults.PhaseCoherence; this.PhaseCoherence = this.CalculatorDefaults.PhaseCoherence;
this.TotalEnergy = this.CalculatorDefaults.TotalEnergy; this.TotalEnergy = this.CalculatorDefaults.TotalEnergy;
this.CutoffType = this.CalculatorDefaults.CutoffType; this.CutoffType = this.CalculatorDefaults.CutoffType;
this.CustomCylindricalCutOffRadius = this.CalculatorDefaults.CustomCylindricalCutOffRadius;
this.CustomCylindricalCutOffHeight = this.CalculatorDefaults.CustomCylindricalCutOffHeight;
end end
end % end %

14
Dipolar-Gas-Simulator/+Simulator/@Calculator/calculateVDk.m
View File

@ -58,6 +58,20 @@ function VDk = calculateVDk(this,Params,Transf,TransfRad,IncludeDDICutOff)
K = sqrt(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2); K = sqrt(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
VDk = (cos(alph).^2-1/3).*(1 + 3*cos(Rcut*K)./(Rcut^2.*K.^2) - 3*sin(Rcut*K)./(Rcut^3.*K.^3)); VDk = (cos(alph).^2-1/3).*(1 + 3*cos(Rcut*K)./(Rcut^2.*K.^2) - 3*sin(Rcut*K)./(Rcut^3.*K.^3));
case 'CustomCylindrical'
alph = acos((Transf.KX*sin(Params.theta)*cos(Params.phi)+Transf.KY*sin(Params.theta)*sin(Params.phi)+Transf.KZ*cos(Params.theta))./sqrt(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2));
alph(1) = pi/2;
VDk = cos(alph).^2-1/3;
% Implementing a custom cylindrical cutoff:
VDr = ifftn(VDk);
VDr = fftshift(VDr);
% Define the cylindrical mask
cylinder_mask = (Transf.X.^2 + Transf.Y.^2 <= this.CustomCylindricalCutOffRadius^2) & (abs(Transf.Z) <= this.CustomCylindricalCutOffHeight / 2);
% Multiply the mask to apply the cut-off
VDr = VDr.*double(cylinder_mask);
VDr = ifftshift(VDr);
VDk = fftn(VDr);
otherwise otherwise
disp('Choose a valid DDI cutoff type!') disp('Choose a valid DDI cutoff type!')
return return

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

@ -55,7 +55,7 @@ classdef DipolarGas < handle & matlab.mixin.Copyable
addParameter(p, 'SimulationMode', 'ImaginaryTimeEvolution',... addParameter(p, 'SimulationMode', 'ImaginaryTimeEvolution',...
@(x) assert(any(strcmpi(x,{'ImaginaryTimeEvolution','RealTimeEvolution'})))); @(x) assert(any(strcmpi(x,{'ImaginaryTimeEvolution','RealTimeEvolution'}))));
addParameter(p, 'CutoffType', 'Cylindrical',... addParameter(p, 'CutoffType', 'Cylindrical',...
@(x) assert(any(strcmpi(x,{'Cylindrical','CylindricalInfiniteZ', 'Spherical'})))); @(x) assert(any(strcmpi(x,{'Cylindrical','CylindricalInfiniteZ', 'Spherical', 'CustomCylindrical'}))));
addParameter(p, 'TimeStepSize', 5E-4,... addParameter(p, 'TimeStepSize', 5E-4,...
@(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); @(x) assert(isnumeric(x) && isscalar(x) && (x > 0)));
addParameter(p, 'MinimumTimeStepSize', 1e-6,... addParameter(p, 'MinimumTimeStepSize', 1e-6,...