From 6cbeb19473b79acd7b11b0a9abeea3a3deaa546c Mon Sep 17 00:00:00 2001 From: Karthik Date: Tue, 1 Apr 2025 01:18:54 +0200 Subject: [PATCH] Backing up of script of CGD implementation. --- Estimations/MinimizeEnergyFunctionalViaCGD.m | 171 +++++++++++++++++++ 1 file changed, 171 insertions(+) create mode 100644 Estimations/MinimizeEnergyFunctionalViaCGD.m diff --git a/Estimations/MinimizeEnergyFunctionalViaCGD.m b/Estimations/MinimizeEnergyFunctionalViaCGD.m new file mode 100644 index 0000000..51d42d6 --- /dev/null +++ b/Estimations/MinimizeEnergyFunctionalViaCGD.m @@ -0,0 +1,171 @@ +function [psi] = minimizeEnergyFunctional(this,psi,Params,Transf,VDk,V,Observ) + + format long; + + % 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 + +%% Helper functions +% Numerical Gradient Calculation using the finite differences method +function J = compute_gradient(psi, Params, Transf, VDk, V) + + % Operators + + % Kinetic energy + KEop = 0.5 * (Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2); + HKin = @(w) real(ifft(KEop.*fft(w))); + + % Trap Potential + HV = @(w) V.*w; + + % Contact interactions + Hint = @(w) (Params.gs*abs(psi).^2).*w; + + % DDIs + frho = fftn(abs(psi).^2); + Phi = real(ifftn(frho.*VDk)); + Hddi = @(w) (Params.gdd*Phi).*w; + + % Quantum fluctuations + Hqf = @(w) (Params.gammaQF*abs(psi).^3).*w; + + H = @(w) HKin(w) + HV(w) + Hint(w) + Hddi(w) + Hqf(w); + + J = H(psi); + +end + +% 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