diff --git a/Dipolar Gas Simulator/+Scripts/test.m b/Dipolar Gas Simulator/+Scripts/test.m index c2679d7..5fa5416 100644 --- a/Dipolar Gas Simulator/+Scripts/test.m +++ b/Dipolar Gas Simulator/+Scripts/test.m @@ -28,11 +28,10 @@ options = Helper.convertstruct2cell(OptionsStruct) clear OptionsStruct sim = Simulator.DipolarGas(options{:}); -calc = Simulator.Calculator(options{:}); pot = Simulator.Potentials(options{:}); %-% Run Simulation %-% -[Params, Transf, psi, V, VDk] = sim.runSimulation(calc); +[Params, Transf, psi, V, VDk] = sim.runSimulation(); %% - Plot numerical grid Plotter.visualizeSpace(Transf) diff --git a/Dipolar Gas Simulator/+Simulator/@Calculator/Calculator.m b/Dipolar Gas Simulator/+Simulator/@Calculator/Calculator.m index 16522d3..93469f0 100644 --- a/Dipolar Gas Simulator/+Simulator/@Calculator/Calculator.m +++ b/Dipolar Gas Simulator/+Simulator/@Calculator/Calculator.m @@ -27,13 +27,23 @@ classdef Calculator < handle & matlab.mixin.Copyable methods function this = Calculator(varargin) - this.ChemicalPotential = this.CalculatorDefaults.ChemicalPotential; - this.EnergyComponents = this.CalculatorDefaults.EnergyComponents; - this.NormalizedResiduals = this.CalculatorDefaults.NormalizedResiduals; - this.OrderParameter = this.CalculatorDefaults.OrderParameter; - this.PhaseCoherence = this.CalculatorDefaults.PhaseCoherence; - this.TotalEnergy = this.CalculatorDefaults.TotalEnergy; - this.CutoffType = this.CalculatorDefaults.CutoffType; + + p = inputParser; + p.KeepUnmatched = true; + addParameter(p, 'CutoffType', this.CalculatorDefaults.CutoffType,... + @(x) any(strcmpi(x,{'Cylindrical','CylindricalInfiniteZ', 'Spherical'}))); + + p.parse(varargin{:}); + + this.ChemicalPotential = this.CalculatorDefaults.ChemicalPotential; + this.EnergyComponents = this.CalculatorDefaults.EnergyComponents; + this.NormalizedResiduals = this.CalculatorDefaults.NormalizedResiduals; + this.OrderParameter = this.CalculatorDefaults.OrderParameter; + this.PhaseCoherence = this.CalculatorDefaults.PhaseCoherence; + this.TotalEnergy = this.CalculatorDefaults.TotalEnergy; + this.CutoffType = p.Results.CutoffType; + this.CalculatorDefaults.CutoffType = this.CutoffType; + end function restoreDefaults(this) diff --git a/Dipolar Gas Simulator/+Simulator/@Calculator/calculateChemicalPotential.m b/Dipolar Gas Simulator/+Simulator/@Calculator/calculateChemicalPotential.m index a7c1c05..96a6a84 100644 --- a/Dipolar Gas Simulator/+Simulator/@Calculator/calculateChemicalPotential.m +++ b/Dipolar Gas Simulator/+Simulator/@Calculator/calculateChemicalPotential.m @@ -1,28 +1,29 @@ -function muchem = calculateChemicalPotential(psi,Params,Transf,VDk,V) +function muchem = calculateChemicalPotential(~,psi,Params,Transf,VDk,V) -%Parameters -normfac = Params.Lx*Params.Ly*Params.Lz/numel(psi); -KEop= 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2); - -% DDIs -frho=fftn(abs(psi).^2); -Phi=real(ifftn(frho.*VDk)); - -Eddi = (Params.gdd*Phi.*abs(psi).^2); - -%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; - -%Contact interactions -Eint = Params.gs*abs(psi).^4; - -%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; \ No newline at end of file + %Parameters + normfac = Params.Lx*Params.Ly*Params.Lz/numel(psi); + KEop= 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2); + + % DDIs + frho=fftn(abs(psi).^2); + Phi=real(ifftn(frho.*VDk)); + + Eddi = (Params.gdd*Phi.*abs(psi).^2); + + %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; + + %Contact interactions + Eint = Params.gs*abs(psi).^4; + + %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; +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 327bfde..92dd32c 100644 --- a/Dipolar Gas Simulator/+Simulator/@Calculator/calculateEnergyComponents.m +++ b/Dipolar Gas Simulator/+Simulator/@Calculator/calculateEnergyComponents.m @@ -1,35 +1,35 @@ -function E = calculateEnergyComponents(psi,Params,Transf,VDk,V) +function E = calculateEnergyComponents(~,psi,Params,Transf,VDk,V) -%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)); - -Eddi = 0.5*Params.gdd*Phi.*abs(psi).^2; -E.Eddi = trapz(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; - -% Potential energy -Epot = V.*abs(psi).^2; -E.Epot = trapz(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; - -%Quantum fluctuations -Eqf = 0.4*Params.gammaQF*abs(psi).^5; -E.Eqf = trapz(Eqf(:))*Transf.dx*Transf.dy*Transf.dz; - -% plot(Transf.x,abs(psi(:,end/2,end/2+1)).^2) \ No newline at end of file + %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)); + + Eddi = 0.5*Params.gdd*Phi.*abs(psi).^2; + E.Eddi = trapz(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; + + % Potential energy + Epot = V.*abs(psi).^2; + E.Epot = trapz(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; + + %Quantum fluctuations + Eqf = 0.4*Params.gammaQF*abs(psi).^5; + E.Eqf = trapz(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 a08e895..134e426 100644 --- a/Dipolar Gas Simulator/+Simulator/@Calculator/calculateNormalizedResiduals.m +++ b/Dipolar Gas Simulator/+Simulator/@Calculator/calculateNormalizedResiduals.m @@ -1,24 +1,25 @@ -function res = calculateNormalizedResiduals(psi,Params,Transf,VDk,V,muchem) +function res = calculateNormalizedResiduals(~,psi,Params,Transf,VDk,V,muchem) -KEop= 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2); - -% DDIs -frho=fftn(abs(psi).^2); -Phi=real(ifftn(frho.*VDk)); - -Eddi = Params.gdd*Phi.*psi; - -%Kinetic energy -Ekin = ifftn(KEop.*fftn(psi)); - -%Potential energy -Epot = V.*psi; - -%Contact interactions -Eint = Params.gs*abs(psi).^2.*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); \ No newline at end of file + KEop= 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2); + + % DDIs + frho=fftn(abs(psi).^2); + Phi=real(ifftn(frho.*VDk)); + + Eddi = Params.gdd*Phi.*psi; + + %Kinetic energy + Ekin = ifftn(KEop.*fftn(psi)); + + %Potential energy + Epot = V.*psi; + + %Contact interactions + Eint = Params.gs*abs(psi).^2.*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); +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 3128b91..b3d96c6 100644 --- a/Dipolar Gas Simulator/+Simulator/@Calculator/calculateNumericalHankelTransform.m +++ b/Dipolar Gas Simulator/+Simulator/@Calculator/calculateNumericalHankelTransform.m @@ -1,39 +1,39 @@ -function VDkSemi = calculateNumericalHankelTransform(this,kr,kz,Rmax,Zmax,Nr) - -% accuracy inputs for numerical integration -if(nargin==5) - Nr = 5e4; -end - -Nz = 64; -farRmultiple = 2000; - -% midpoint grids for the integration over 0 0))); addParameter(p, 'SimulationMode', 'ImaginaryTimeEvolution',... @(x) any(strcmpi(x,{'ImaginaryTimeEvolution','RealTimeEvolution'}))); + addParameter(p, 'CutoffType', 'Cylindrical',... + @(x) any(strcmpi(x,{'Cylindrical','CylindricalInfiniteZ', 'Spherical'}))); addParameter(p, 'TimeStep', 5E-4,... @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); addParameter(p, 'SimulationTime', 2e6,... @@ -81,6 +85,8 @@ classdef DipolarGas < handle & matlab.mixin.Copyable this.DoSave = p.Results.SaveData; this.SaveDirectory = p.Results.SaveDirectory; + this.Calculator = Simulator.Calculator('CutoffType', p.Results.CutoffType); + switch this.SimulationMode case "ImaginaryTimeEvolution" % Development In progress diff --git a/Dipolar Gas Simulator/+Simulator/@DipolarGas/initialize.m b/Dipolar Gas Simulator/+Simulator/@DipolarGas/initialize.m index d7d5eb0..ef860d5 100644 --- a/Dipolar Gas Simulator/+Simulator/@DipolarGas/initialize.m +++ b/Dipolar Gas Simulator/+Simulator/@DipolarGas/initialize.m @@ -1,22 +1,20 @@ -function [psi,V,VDk] = initialize(this,calcObj,Params,Transf,TransfRad) - -format long -X = Transf.X; Y = Transf.Y; Z = Transf.Z; - -% == Trap potential == % -V = 0.5*(Params.gx.*X.^2+Params.gy.*Y.^2+Params.gz*Z.^2); - -% == Calculating the DDIs == % -if isfile(strcat(this.SaveDirectory, '/VDk_M.mat')) - VDk = load(sprintf(strcat(this.SaveDirectory, '/VDk_M.mat'))); - VDk = VDk.VDk; -else - VDk = calcObj.calculateVDCutoff(Params,Transf,TransfRad); - save(sprintf(strcat(this.SaveDirectory, '/VDk_M.mat')),'VDk'); -end -fprintf('Computed and saved DDI potential in Fourier space with %s cutoff.', calcObj.CutoffType) - -% == Setting up the initial wavefunction == % -psi = this.setupWavefunction(Params,Transf); - +function [psi,V,VDk] = initialize(this,Params,Transf,TransfRad) + format long + X = Transf.X; Y = Transf.Y; Z = Transf.Z; + + % == Trap potential == % + V = 0.5*(Params.gx.*X.^2+Params.gy.*Y.^2+Params.gz*Z.^2); + + % == Calculating the DDIs == % + if isfile(strcat(this.SaveDirectory, '/VDk_M.mat')) + VDk = load(sprintf(strcat(this.SaveDirectory, '/VDk_M.mat'))); + VDk = VDk.VDk; + else + VDk = this.Calculator.calculateVDCutoff(Params,Transf,TransfRad); + save(sprintf(strcat(this.SaveDirectory, '/VDk_M.mat')),'VDk'); + end + fprintf('Computed and saved DDI potential in Fourier space with %s cutoff.', this.Calculator.CutoffType) + + % == Setting up the initial wavefunction == % + psi = this.setupWavefunction(Params,Transf); end \ No newline at end of file diff --git a/Dipolar Gas Simulator/+Simulator/@DipolarGas/runSimulation.m b/Dipolar Gas Simulator/+Simulator/@DipolarGas/runSimulation.m index 74ef82b..aba48d2 100644 --- a/Dipolar Gas Simulator/+Simulator/@DipolarGas/runSimulation.m +++ b/Dipolar Gas Simulator/+Simulator/@DipolarGas/runSimulation.m @@ -1,4 +1,4 @@ -function [Params, Transf, psi,V,VDk] = runSimulation(this,calcObj) +function [Params, Transf, psi,V,VDk] = runSimulation(this) % --- Obtain simulation parameters --- [Params] = this.setupParameters(); @@ -9,15 +9,15 @@ function [Params, Transf, psi,V,VDk] = runSimulation(this,calcObj) % --- Initialize --- mkdir(sprintf(this.SaveDirectory)) - [psi,V,VDk] = this.initialize(calcObj,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; % --- Job Settings --- - % njob = 6; - % mkdir(sprintf('./Data/Run_%03i',njob)) + njob = 6; + mkdir(sprintf('./Data/Run_%03i',njob)) % --- Run Simulation --- % [psi] = this.solver(psi,Params,Transf,VDk,V,njob,t_idx,Observ); diff --git a/Dipolar Gas Simulator/+Simulator/@DipolarGas/setupParameters.m b/Dipolar Gas Simulator/+Simulator/@DipolarGas/setupParameters.m index bf6114e..e0f516d 100644 --- a/Dipolar Gas Simulator/+Simulator/@DipolarGas/setupParameters.m +++ b/Dipolar Gas Simulator/+Simulator/@DipolarGas/setupParameters.m @@ -1,92 +1,91 @@ function [Params] = setupParameters(this) - -CONSTANTS = Helper.PhysicsConstants; -hbar = CONSTANTS.PlanckConstantReduced; % [J.s] -kbol = CONSTANTS.BoltzmannConstant; % [J/K] -mu0 = CONSTANTS.VacuumPermeability; % [N/A^2] -muB = CONSTANTS.BohrMagneton; % [J/T] -a0 = CONSTANTS.BohrRadius; % [m] -m0 = CONSTANTS.AtomicMassUnit; % [kg] -w0 = 2*pi*100; % 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 -Params.Nx = this.NumberOfGridPoints(1); -Params.Ny = this.NumberOfGridPoints(2); -Params.Nz = this.NumberOfGridPoints(3); - -% Dimensions (in units of l0) -Params.Lx = this.Dimensions(1); -Params.Ly = this.Dimensions(2); -Params.Lz = this.Dimensions(3); - -% Masses -Params.m = CONSTANTS.Dy164Mass; -l0 = sqrt(hbar/(Params.m*w0)); % Defining a harmonic oscillator length - -% Atom numbers -Params.N = this.NumberOfAtoms; - -% Dipole angle -Params.theta = this.DipolarPolarAngle; % pi/2 dipoles along x, theta=0 dipoles along z -Params.phi = this.DipolarAzimuthAngle; - -% Dipole lengths (units of muB) -Params.mu = CONSTANTS.DyMagneticMoment; - -% Scattering lengths -Params.as = this.ScatteringLength*a0; - -% Trapping frequencies -Params.wx = 2*pi*this.TrapFrequencies(1); -Params.wy = 2*pi*this.TrapFrequencies(2); -Params.wz = 2*pi*this.TrapFrequencies(3); - -% Stochastic GPE -Params.gamma_S = 7.5*10^(-3); % gamma for the stochastic GPE -Params.muchem = 12.64*Params.wz/w0; % fixing the chemical potential for the stochastic GPE - -Params.Etol = this.EnergyTolerance; % Tolerances -Params.cut_off = this.SimulationTime; % sometimes the imaginary time gets a little stuck - % even though the solution is good, this just stops it going on forever -Params.mindt = this.MinimumTimeStep; % Minimum size for a time step using adaptive dt - -% ================ Parameters defined by those above ================ % - -% Contact interaction strength (units of l0/m) -Params.gs = 4*pi*Params.as/l0; - -% Dipole lengths -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 - -% Trap gamma -Params.gx = (Params.wx/w0)^2; -Params.gy = (Params.wy/w0)^2; -Params.gz = (Params.wz/w0)^2; - -% == Calculate LHY correction == % -eps_dd = Params.add/Params.as; -if eps_dd == 0 - Q5 = 1; -elseif eps_dd == 1 - Q5 = 3*sqrt(3)/2; -else - yeps = (1-eps_dd)/(3*eps_dd); - Q5 = (3*eps_dd)^(5/2)*( (8+26*yeps+33*yeps^2)*sqrt(1+yeps) + 15*yeps^3*log((1+sqrt(1+yeps))/sqrt(yeps)) )/48; - Q5 = real(Q5); -end - -Params.gammaQF = 128/3*sqrt(pi*(Params.as/l0)^5)*Q5; - -% Loading the rest into Params -Params.hbar = hbar; -Params.kbol = kbol; -Params.mu0 = mu0; -Params.muB = muB; -Params.a0 = a0; -Params.w0 = w0; -Params.l0 = l0; + CONSTANTS = Helper.PhysicsConstants; + hbar = CONSTANTS.PlanckConstantReduced; % [J.s] + kbol = CONSTANTS.BoltzmannConstant; % [J/K] + mu0 = CONSTANTS.VacuumPermeability; % [N/A^2] + muB = CONSTANTS.BohrMagneton; % [J/T] + a0 = CONSTANTS.BohrRadius; % [m] + m0 = CONSTANTS.AtomicMassUnit; % [kg] + w0 = 2*pi*100; % 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 + Params.Nx = this.NumberOfGridPoints(1); + Params.Ny = this.NumberOfGridPoints(2); + Params.Nz = this.NumberOfGridPoints(3); + + % Dimensions (in units of l0) + Params.Lx = this.Dimensions(1); + Params.Ly = this.Dimensions(2); + Params.Lz = this.Dimensions(3); + + % Masses + Params.m = CONSTANTS.Dy164Mass; + l0 = sqrt(hbar/(Params.m*w0)); % Defining a harmonic oscillator length + + % Atom numbers + Params.N = this.NumberOfAtoms; + + % Dipole angle + Params.theta = this.DipolarPolarAngle; % pi/2 dipoles along x, theta=0 dipoles along z + Params.phi = this.DipolarAzimuthAngle; + + % Dipole lengths (units of muB) + Params.mu = CONSTANTS.DyMagneticMoment; + + % Scattering lengths + Params.as = this.ScatteringLength*a0; + + % Trapping frequencies + Params.wx = 2*pi*this.TrapFrequencies(1); + Params.wy = 2*pi*this.TrapFrequencies(2); + Params.wz = 2*pi*this.TrapFrequencies(3); + + % Stochastic GPE + Params.gamma_S = 7.5*10^(-3); % gamma for the stochastic GPE + Params.muchem = 12.64*Params.wz/w0; % fixing the chemical potential for the stochastic GPE + + Params.Etol = this.EnergyTolerance; % Tolerances + Params.cut_off = this.SimulationTime; % sometimes the imaginary time gets a little stuck + % even though the solution is good, this just stops it going on forever + Params.mindt = this.MinimumTimeStep; % Minimum size for a time step using adaptive dt + + % ================ Parameters defined by those above ================ % + + % Contact interaction strength (units of l0/m) + Params.gs = 4*pi*Params.as/l0; + + % Dipole lengths + 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 + + % Trap gamma + Params.gx = (Params.wx/w0)^2; + Params.gy = (Params.wy/w0)^2; + Params.gz = (Params.wz/w0)^2; + + % == Calculate LHY correction == % + eps_dd = Params.add/Params.as; + if eps_dd == 0 + Q5 = 1; + elseif eps_dd == 1 + Q5 = 3*sqrt(3)/2; + else + yeps = (1-eps_dd)/(3*eps_dd); + Q5 = (3*eps_dd)^(5/2)*( (8+26*yeps+33*yeps^2)*sqrt(1+yeps) + 15*yeps^3*log((1+sqrt(1+yeps))/sqrt(yeps)) )/48; + Q5 = real(Q5); + end + + Params.gammaQF = 128/3*sqrt(pi*(Params.as/l0)^5)*Q5; + + % Loading the rest into Params + Params.hbar = hbar; + Params.kbol = kbol; + Params.mu0 = mu0; + Params.muB = muB; + Params.a0 = a0; + Params.w0 = w0; + Params.l0 = l0; end \ No newline at end of file diff --git a/Dipolar Gas Simulator/+Simulator/@DipolarGas/setupSpace.m b/Dipolar Gas Simulator/+Simulator/@DipolarGas/setupSpace.m index 56f7f7c..56996bd 100644 --- a/Dipolar Gas Simulator/+Simulator/@DipolarGas/setupSpace.m +++ b/Dipolar Gas Simulator/+Simulator/@DipolarGas/setupSpace.m @@ -1,33 +1,33 @@ function [Transf] = setupSpace(this,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; - -% 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); - -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); - -[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.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); + + 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); + + [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/setupSpaceRadial.m b/Dipolar Gas Simulator/+Simulator/@DipolarGas/setupSpaceRadial.m index f9349e2..2686b0e 100644 --- a/Dipolar Gas Simulator/+Simulator/@DipolarGas/setupSpaceRadial.m +++ b/Dipolar Gas Simulator/+Simulator/@DipolarGas/setupSpaceRadial.m @@ -1,51 +1,50 @@ -function [Transf] = setupSpaceRadial(this,Params,morder) - -Params.Lr = 0.5*min(Params.Lx,Params.Ly); -Params.Nr = max(Params.Nx,Params.Ny); - -Zmax = 0.5*Params.Lz; -Rmax = Params.Lr; -Nz = Params.Nz; -Nr = Params.Nr; - -if(nargin==2) - morder=0; %only do Bessel J0 -end - -% Fourier grids -z=linspace(-Zmax,Zmax,Nz+1); -z=z(1:end-1); -dz=z(2)-z(1); -Kmax=Nz*2*pi/(4*Zmax); -kz=linspace(-Kmax,Kmax,Nz+1); -kz=kz(1:end-1); - -% Hankel grids and transform -H = hankelmatrix(morder,Rmax,Nr); -r=H.r(:); -kr=H.kr(:); -T = diag(H.J/H.kmax)*H.T*diag(Rmax./H.J)*dz*(2*pi); -Tinv = diag(H.J./Rmax)*H.T'*diag(H.kmax./H.J)/dz/(2*pi); -wr=H.wr; -wk=H.wk; -% H.T'*diag(H.J/H.vmax)*H.T*diag(Rmax./H.J) - -[Transf.R,Transf.Z]=ndgrid(r,z); -[Transf.KR,Transf.KZ]=ndgrid(kr,kz); -Transf.T=T; -Transf.Tinv=Tinv; -Transf.r=r; -Transf.kr=kr; -Transf.z=z; -Transf.kz=kz; -Transf.wr=wr; -Transf.wk=wk; -Transf.Rmax=Rmax; -Transf.Zmax=Zmax; -Transf.dz=z(2)-z(1); -Transf.dkz=kz(2)-kz(1); -%b1=Transf; +function [Transf] = setupSpaceRadial(~,Params,morder) + Params.Lr = 0.5*min(Params.Lx,Params.Ly); + Params.Nr = max(Params.Nx,Params.Ny); + + Zmax = 0.5*Params.Lz; + Rmax = Params.Lr; + Nz = Params.Nz; + Nr = Params.Nr; + + if(nargin==2) + morder=0; %only do Bessel J0 + end + + % Fourier grids + z=linspace(-Zmax,Zmax,Nz+1); + z=z(1:end-1); + dz=z(2)-z(1); + Kmax=Nz*2*pi/(4*Zmax); + kz=linspace(-Kmax,Kmax,Nz+1); + kz=kz(1:end-1); + + % Hankel grids and transform + H = hankelmatrix(morder,Rmax,Nr); + r=H.r(:); + kr=H.kr(:); + T = diag(H.J/H.kmax)*H.T*diag(Rmax./H.J)*dz*(2*pi); + Tinv = diag(H.J./Rmax)*H.T'*diag(H.kmax./H.J)/dz/(2*pi); + wr=H.wr; + wk=H.wk; + % H.T'*diag(H.J/H.vmax)*H.T*diag(Rmax./H.J) + + [Transf.R,Transf.Z]=ndgrid(r,z); + [Transf.KR,Transf.KZ]=ndgrid(kr,kz); + Transf.T=T; + Transf.Tinv=Tinv; + Transf.r=r; + Transf.kr=kr; + Transf.z=z; + Transf.kz=kz; + Transf.wr=wr; + Transf.wk=wk; + Transf.Rmax=Rmax; + Transf.Zmax=Zmax; + Transf.dz=z(2)-z(1); + Transf.dkz=kz(2)-kz(1); + function s_HT = hankelmatrix(order,rmax,Nr,eps_roots) %HANKEL_MATRIX: Generates data to use for Hankel Transforms % @@ -103,9 +102,9 @@ function s_HT = hankelmatrix(order,rmax,Nr,eps_roots) % See also bessel_zeros, besselj if (~exist('eps_roots', 'var')||isemtpy(eps_roots)) - s_HT.eps_roots = 1e-6; + s_HT.eps_roots = 1e-6; else - s_HT.eps_roots = eps_roots; + s_HT.eps_roots = eps_roots; end s_HT.order = order; @@ -166,146 +165,146 @@ mu3 = mu^3; mu4 = mu^4; if (d<3) - p = 7*mu - 31; - p0 = mu - 1; - if ((1+p)==p) - p1 = 0; - q1 = 0; - else - p1 = 4*(253*mu2 - 3722*mu+17869)*p0/(15*p); - q1 = 1.6*(83*mu2 - 982*mu + 3779)/p; - end + p = 7*mu - 31; + p0 = mu - 1; + if ((1+p)==p) + p1 = 0; + q1 = 0; + else + p1 = 4*(253*mu2 - 3722*mu+17869)*p0/(15*p); + q1 = 1.6*(83*mu2 - 982*mu + 3779)/p; + end else - p = 7*mu2 + 82*mu - 9; - p0 = mu + 3; - if ((p+1)==1) - p1 = 0; - q1 = 0; - else - p1 = (4048*mu4 + 131264*mu3 - 221984*mu2 - 417600*mu + 1012176)/(60*p); - q1 = 1.6*(83*mu3 + 2075*mu2 - 3039*mu + 3537)/p; - end + p = 7*mu2 + 82*mu - 9; + p0 = mu + 3; + if ((p+1)==1) + p1 = 0; + q1 = 0; + else + p1 = (4048*mu4 + 131264*mu3 - 221984*mu2 - 417600*mu + 1012176)/(60*p); + q1 = 1.6*(83*mu3 + 2075*mu2 - 3039*mu + 3537)/p; + end end if (d==1)|(d==4) - t = .25; + t = .25; else - t = .75; + t = .75; end tt = 4*t; if (d<3) - pp1 = 5/48; - qq1 = -5/36; + pp1 = 5/48; + qq1 = -5/36; else - pp1 = -7/48; - qq1 = 35/288; + pp1 = -7/48; + qq1 = 35/288; end y = .375*pi; if (a>=3) - bb = a^(-2/3); + bb = a^(-2/3); else - bb = 1; + bb = 1; end a1 = 3*a - 8; % psi = (.5*a + .25)*pi; for s=1:n - if ((a==0)&(s==1)&(d==3)) - x = 0; - j = 0; - else - if (s>=a1) - b = (s + .5*a - t)*pi; - c = .015625/(b^2); - x = b - .125*(p0 - p1*c)/(b*(1 - q1*c)); - else - if (s==1) - switch (d) - case (1) - x = -2.33811; - case (2) - x = -1.17371; - case (3) - x = -1.01879; - otherwise - x = -2.29444; - end - else - x = y*(4*s - tt); - v = x^(-2); - x = -x^(2/3) * (1 + v*(pp1 + qq1*v)); - end - u = x*bb; - v = fi(2/3 * (-u)^1.5); - w = 1/cos(v); - xx = 1 - w^2; - c = sqrt(u/xx); - if (d<3) - x = w*(a + c*(-5/u - c*(6 - 10/xx))/(48*a*u)); - else - x = w*(a + c*(7/u + c*(18 - 14/xx))/(48*a*u)); - end - end - j = 0; - + if ((a==0)&(s==1)&(d==3)) + x = 0; + j = 0; + else + if (s>=a1) + b = (s + .5*a - t)*pi; + c = .015625/(b^2); + x = b - .125*(p0 - p1*c)/(b*(1 - q1*c)); + else + if (s==1) + switch (d) + case (1) + x = -2.33811; + case (2) + x = -1.17371; + case (3) + x = -1.01879; + otherwise + x = -2.29444; + end + else + x = y*(4*s - tt); + v = x^(-2); + x = -x^(2/3) * (1 + v*(pp1 + qq1*v)); + end + u = x*bb; + v = fi(2/3 * (-u)^1.5); + w = 1/cos(v); + xx = 1 - w^2; + c = sqrt(u/xx); + if (d<3) + x = w*(a + c*(-5/u - c*(6 - 10/xx))/(48*a*u)); + else + x = w*(a + c*(7/u + c*(18 - 14/xx))/(48*a*u)); + end + end + j = 0; + while ((j==0)|((j<5)&(abs(w/x)>e))) - xx = x^2; - x4 = x^4; - a2 = aa - xx; - r0 = bessr(d, a, x); - j = j+1; - if (d<3) - u = r0; - w = 6*x*(2*a + 1); - p = (1 - 4*a2)/w; - q = (4*(xx-mu) - 2 - 12*a)/w; - else - u = -xx*r0/a2; - v = 2*x*a2/(3*(aa+xx)); - w = 64*a2^3; - q = 2*v*(1 + mu2 + 32*mu*xx + 48*x4)/w; - p = v*(1 + (40*mu*xx + 48*x4 - mu2)/w); - end - w = u*(1 + p*r0)/(1 + q*r0); - x = x+w; - end - z(s) = x; - end + xx = x^2; + x4 = x^4; + a2 = aa - xx; + r0 = bessr(d, a, x); + j = j+1; + if (d<3) + u = r0; + w = 6*x*(2*a + 1); + p = (1 - 4*a2)/w; + q = (4*(xx-mu) - 2 - 12*a)/w; + else + u = -xx*r0/a2; + v = 2*x*a2/(3*(aa+xx)); + w = 64*a2^3; + q = 2*v*(1 + mu2 + 32*mu*xx + 48*x4)/w; + p = v*(1 + (40*mu*xx + 48*x4 - mu2)/w); + end + w = u*(1 + p*r0)/(1 + q*r0); + x = x+w; + end + z(s) = x; + end end function FI = fi(y) - c1 = 1.570796; - if (~y) - FI = 0; - elseif (y>1e5) - FI = c1; - else - if (y<1) - p = (3*y)^(1/3); - pp = p^2; - p = p*(1 + pp*(pp*(27 - 2*pp) - 210)/1575); - else - p = 1/(y + c1); - pp = p^2; - p = c1 - p*(1 + pp*(2310 + pp*(3003 + pp*(4818 + pp*(8591 + pp*16328))))/3465); - end - pp = (y+p)^2; - r = (p - atan(p+y))/pp; - FI = p - (1+pp)*r*(1 + r/(p+y)); - end + c1 = 1.570796; + if (~y) + FI = 0; + elseif (y>1e5) + FI = c1; + else + if (y<1) + p = (3*y)^(1/3); + pp = p^2; + p = p*(1 + pp*(pp*(27 - 2*pp) - 210)/1575); + else + p = 1/(y + c1); + pp = p^2; + p = c1 - p*(1 + pp*(2310 + pp*(3003 + pp*(4818 + pp*(8591 + pp*16328))))/3465); + end + pp = (y+p)^2; + r = (p - atan(p+y))/pp; + FI = p - (1+pp)*r*(1 + r/(p+y)); + end return function Jr = bessr(d,a,x) - switch (d) - case (1) - Jr = besselj(a, x)./besselj(a+1, x); - case (2) - Jr = bessely(a, x)./bessely(a+1, x); - case (3) - Jr = a./x - besselj(a+1, x)./besselj(a, x); - otherwise - Jr = a./x - bessely(a+1, x)./bessely(a, x); - end + switch (d) + case (1) + Jr = besselj(a, x)./besselj(a+1, x); + case (2) + Jr = bessely(a, x)./bessely(a+1, x); + case (3) + Jr = a./x - besselj(a+1, x)./besselj(a, x); + otherwise + Jr = a./x - bessely(a+1, x)./bessely(a, x); + end return \ 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 e16cf48..f4c73be 100644 --- a/Dipolar Gas Simulator/+Simulator/@DipolarGas/setupWavefunction.m +++ b/Dipolar Gas Simulator/+Simulator/@DipolarGas/setupWavefunction.m @@ -1,25 +1,25 @@ -function [psi] = setupWavefunction(this,Params,Transf) +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; - -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; -psi = psi/sqrt(cur_norm); - -% add some noise -r = normrnd(0,1,size(X)); -theta = rand(size(X)); -noise = r.*exp(2*pi*1i*theta); -psi = psi + 0.01*noise; - -Norm = trapz(abs(psi(:)).^2)*Transf.dx*Transf.dy*Transf.dz; -psi = sqrt(Params.N)*psi/sqrt(Norm); + 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 = 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; + psi = psi/sqrt(cur_norm); + + % add some noise + r = normrnd(0,1,size(X)); + theta = rand(size(X)); + noise = r.*exp(2*pi*1i*theta); + psi = psi + 0.01*noise; + + Norm = trapz(abs(psi(:)).^2)*Transf.dx*Transf.dy*Transf.dz; + psi = sqrt(Params.N)*psi/sqrt(Norm); end \ No newline at end of file diff --git a/Dipolar Gas Simulator/+Simulator/@DipolarGas/solver.m b/Dipolar Gas Simulator/+Simulator/@DipolarGas/solver.m index bc02386..088f6b4 100644 --- a/Dipolar Gas Simulator/+Simulator/@DipolarGas/solver.m +++ b/Dipolar Gas Simulator/+Simulator/@DipolarGas/solver.m @@ -1,90 +1,89 @@ -function [psi] = solver(this,psi,Params,Transf,VDk,V,njob,t_idx,Observ) - -set(0,'defaulttextInterpreter','latex') -set(groot, 'defaultAxesTickLabelInterpreter','latex'); set(groot, 'defaultLegendInterpreter','latex'); - -dt=-1j*abs(this.TimeStep); - -KEop= 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2); -Observ.residual = 1; Observ.res = 1; - -muchem = Simulator.ChemicalPotential(psi,Params,Transf,VDk,V); -AdaptIdx = 0; - -while t_idx < Params.cut_off - %kin - psi = fftn(psi); - psi = psi.*exp(-0.5*1i*dt*KEop); - psi = ifftn(psi); +function [psi] = propagate(this,psi,Params,Transf,VDk,V,njob,t_idx,Observ) + set(0,'defaulttextInterpreter','latex') + set(groot, 'defaultAxesTickLabelInterpreter','latex'); set(groot, 'defaultLegendInterpreter','latex'); - %DDI - frho = fftn(abs(psi).^2); - Phi = real(ifftn(frho.*VDk)); + dt=-1j*abs(this.TimeStep); - %Real-space - psi = psi.*exp(-1i*dt*(V + Params.gs*abs(psi).^2 + Params.gammaQF*abs(psi).^3 + Params.gdd*Phi - muchem)); - - %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); - - muchem = Simulator.ChemicalPotential(psi,Params,Transf,VDk,V); - - if mod(t_idx,1000) == 0 - - %Change in Energy - E = Simulator.TotalEnergy(psi,Params,Transf,VDk,V); - E = E/Norm; - Observ.EVec = [Observ.EVec E]; - - %Chemical potential - Observ.mucVec = [Observ.mucVec muchem]; - - %Normalized residuals - res = Simulator.NormalizedResiduals(psi,Params,Transf,VDk,V,muchem); - Observ.residual = [Observ.residual res]; + KEop= 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2); + Observ.residual = 1; Observ.res = 1; + + muchem = this.Calculator.ChemicalPotential(psi,Params,Transf,VDk,V); + AdaptIdx = 0; + + while t_idx < Params.cut_off + %kin + psi = fftn(psi); + psi = psi.*exp(-0.5*1i*dt*KEop); + psi = ifftn(psi); - Observ.res_idx = Observ.res_idx + 1; - - save(sprintf('./Data/Run_%03i/psi_gs.mat',njob),'psi','muchem','Observ','t_idx','Transf','Params','VDk','V'); + %DDI + frho = fftn(abs(psi).^2); + Phi = real(ifftn(frho.*VDk)); - %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 AdaptIdx > 4 && abs(dt) > Params.mindt - dt = dt / 2; - fprintf('Time step changed to '); disp(dt); - AdaptIdx = 0; - elseif AdaptIdx > 4 && abs(dt) < Params.mindt - break + %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); + + %Renorm + Norm = trapz(abs(psi(:)).^2)*Transf.dx*Transf.dy*Transf.dz; + psi = sqrt(Params.N)*psi/sqrt(Norm); + + muchem = this.Calculator.ChemicalPotential(psi,Params,Transf,VDk,V); + + if mod(t_idx,1000) == 0 + + %Change in Energy + E = this.Calculator.TotalEnergy(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.NormalizedResiduals(psi,Params,Transf,VDk,V,muchem); + Observ.residual = [Observ.residual res]; + + Observ.res_idx = Observ.res_idx + 1; + + save(sprintf('./Data/Run_%03i/psi_gs.mat',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 AdaptIdx > 4 && abs(dt) > Params.mindt + dt = dt / 2; + fprintf('Time step changed to '); disp(dt); + AdaptIdx = 0; + elseif AdaptIdx > 4 && abs(dt) < Params.mindt + break + else + AdaptIdx = AdaptIdx + 1; + end else - AdaptIdx = AdaptIdx + 1; + AdaptIdx = 0; end - else - AdaptIdx = 0; end + if any(isnan(psi(:))) + disp('NaNs encountered!') + break + end + t_idx=t_idx+1; end - if any(isnan(psi(:))) - disp('NaNs encountered!') - break - end - t_idx=t_idx+1; -end - -%Change in Energy -E = Simulator.TotalEnergy(psi,Params,Transf,VDk,V); -E = E/Norm; -Observ.EVec = [Observ.EVec E]; - -% Phase coherence -[PhaseC] = Simulator.PhaseCoherence(psi,Transf,Params); -Observ.PCVec = [Observ.PCVec PhaseC]; - -Observ.res_idx = Observ.res_idx + 1; -save(sprintf('./Data/Run_%03i/psi_gs.mat',njob),'psi','muchem','Observ','t_idx','Transf','Params','VDk','V'); + + %Change in Energy + E = this.Calculator.TotalEnergy(psi,Params,Transf,VDk,V); + E = E/Norm; + Observ.EVec = [Observ.EVec E]; + + % Phase coherence + [PhaseC] = this.Calculator.PhaseCoherence(psi,Transf,Params); + Observ.PCVec = [Observ.PCVec PhaseC]; + + Observ.res_idx = Observ.res_idx + 1; + save(sprintf('./Data/Run_%03i/psi_gs.mat',njob),'psi','muchem','Observ','t_idx','Transf','Params','VDk','V'); end \ No newline at end of file