SelfConsistent.m

Go to the documentation of this file.

%%    Eotvos Quantum Transport Utilities - SelfConsistent
%    Copyright (C) 2009-2017 Peter Rakyta, Ph.D.
%
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program.  If not, see http://www.gnu.org/licenses/.
%
%> @addtogroup basic Basic Functionalities
%> @{
%> @file SelfConsistent.m
%> @brief A class to evaluate self-consistent potentials depending on electron (spin resolved) densities.
%> @} 
%> @brief A class to evaluate self-consistent potentials depending on electron (spin resolved) densities.
%> @Available
%> EQuUs v4.9 or later
%%
classdef SelfConsistent < handle & Messages

    properties ( Access = public )
%> An array containing the current density
        density
%> structure describing the geometry of the sites
        junction_sites
%> The calculated filling factor
        filling_factor 
%> The prescribed magnezitaion direction (+1/-1)
        magnetizaion_direction
%> The current magnetization
        magnetization        
%> Logical value. Set true to export graphical results for each iteration %NEW
        plotIterations
%> The function handle for which the self-consistent iterations are performed: y = fhandle( qvec ), where y is an MxN matrix, and qvec is a vector of transverse momentums containing of M elements.      
        fhandle
    end
    
    
    properties ( Access = protected )
%> false to supress output messages
        silent
 %> Tolerance in the self-consistent iterations
        tolerance
%> Maximal number of iterations
        MaxIter
%>  Logical value. If true, the iteration cintinues with the loaded data, or false (default) otherwise.
        resume   
%>  The filename to export the calculated data.
        outFileName  
%>  A list containing previously calculated densities
        densities_iterations
%>  A list containing the difference between previously calculated densities
        delta_densities_iterations        
%> A list containing previously calculated magnetizations
        magnetization_iterations
%> Scalar value. The weith of the newly calculated density in the iterations.
        newIterationWeight        
%> list of optional parameters (see http://www.mathworks.com/help/matlab/ref/varargin.html for details)
        varargin         
    end    



methods (Access=public)
    
%% Contructor of the class
%> @brief Constructor of the class.
%> @param Opt An instance of the structure #Opt.
%> @param interpq List of transverse momentum points at which #fhandle is calculated via interpolation.
%> @param varargin Cell array of optional parameters. For details see #InputParsing.
%> @return An instance of the class
    function obj = SelfConsistent( Opt,  varargin )        

        obj = obj@Messages( Opt );
        obj.varargin = varargin;

        obj.Initialize();       
        
        
        
    end
    
    
%% runSelfConsistentIterations
%> @brief Runs the self-consistent iterations to obtain the elemets of the self-consistent densities for spin-up and spin-down electrons
%> @param fhandle The function handle to be calculated over the transverse momentum points: y = fhandle( qvec ), where y is an MxN matrix, and qvec is a vector of transverse momentums containing of M elements.
    function runSelfConsistentIterations( obj, fhandle )
        
        obj.display( '*************************************')
        obj.display( 'Starting self-consistent iterations.', 1 );
        obj.display( '*************************************')
        
        if obj.resume && ~isempty(obj.outFileName)
            
            try      
                HandleForLoad = LoadFromFile( obj.outFileName );
                loaded_data = HandleForLoad.LoadVariable('obj', 'NoEmpty', 'On');           
                HandleForLoad.ClearLoadedVariables();
                
                meta_data = metaclass(loaded_data);
        
                for idx = 1:length(meta_data.PropertyList)
                    prop_name = meta_data.PropertyList(idx).Name;
                    obj.(prop_name) = loaded_data.(prop_name);
                end         
        
            catch err
               obj.display( 'EQuUs:SelfConsitent:runSelfConsistentIterations: Failed to load data. Restarting calculations.', 1) 
            end
        end
        
        %> run the first iteration
        if isempty( obj.density )
            obj.Initialize();
            [obj.density, obj.junction_sites] = feval(fhandle);
        end
        

        %> do the rest of the iterations        
        obj.filling_factor = sum(obj.density)/numel(obj.density);
        for idx = 1:obj.MaxIter
            % break the iterations if convergence is reached
            if obj.ConvergenceTest()
                break;
            end
            
            % densities for the next iterations
            obj.NewIterationInput();          
            
            % calculate the current densities
            [obj.density, obj.junction_sites] = feval(fhandle);
            
%             scatter_density = obj.density(obj.junction_sites.Scatter.site_indexes);
%             figure
%             hold on
%             plot( obj.junction_sites.Scatter.coordinates.x(obj.junction_sites.Scatter.coordinates.spinup), scatter_density(obj.junction_sites.Scatter.coordinates.spinup))
%             plot( obj.junction_sites.Scatter.coordinates.x(~obj.junction_sites.Scatter.coordinates.spinup), scatter_density(~obj.junction_sites.Scatter.coordinates.spinup), 'r')            
            
            obj.filling_factor   =  sum(obj.density)/numel(obj.density);
            
            %> reverting the magnetization direction if necessary
            if obj.Opt.Spin
                scatter_density = obj.density(obj.junction_sites.Scatter.site_indexes);
                obj.magnetization = sum(scatter_density(obj.junction_sites.Scatter.coordinates.spinup) - scatter_density(~obj.junction_sites.Scatter.coordinates.spinup));
                if ~isempty(obj.magnetizaion_direction)
                    if sign(obj.magnetization) ~= sign(obj.magnetizaion_direction)
                        density_tmp = scatter_density;
                        scatter_density(obj.junction_sites.Scatter.coordinates.spinup) = density_tmp(~obj.junction_sites.Scatter.coordinates.spinup);
                        scatter_density(~obj.junction_sites.Scatter.coordinates.spinup) = density_tmp(obj.junction_sites.Scatter.coordinates.spinup);
                        obj.density(obj.junction_sites.Scatter.site_indexes) = scatter_density;
                        obj.magnetization = -obj.magnetization;
                    end
                end
            end            
            
            %> calculating the difference compared to the privous density
            if ~isempty( obj.densities_iterations{1} )
                obj.delta_densities_iterations(1) = norm(obj.densities_iterations{1}-obj.density)/norm(obj.densities_iterations{1});
            end
            
            
            %> creating and displaying output messaage
            message = ['Iteration = ', num2str(idx), ', filling factor = ', num2str(obj.filling_factor), ', density delta = ', num2str(obj.delta_densities_iterations(1))];
            if obj.Opt.Spin
               message = [message, ', Magnetization = ', num2str(obj.magnetization_iterations(1)), ', magnetization delta = ', num2str(abs((obj.magnetization_iterations(1) - obj.magnetization)/obj.magnetization))];
            end
            obj.display(message, 1);
            
            if obj.plotIterations
                scatter_density = obj.density(obj.junction_sites.Scatter.site_indexes);
                figure1 = figure;
                hold on
                plot( obj.junction_sites.Scatter.coordinates.x(obj.junction_sites.Scatter.coordinates.spinup), scatter_density(obj.junction_sites.Scatter.coordinates.spinup), 'b+')
                plot( obj.junction_sites.Scatter.coordinates.x(~obj.junction_sites.Scatter.coordinates.spinup), scatter_density(~obj.junction_sites.Scatter.coordinates.spinup), 'ro')
                xlabel('x')
                ylabel('Density')
                print('-depsc2', ['Iteration_', num2str(idx), '.eps'])
                close(figure1);
            end
            
            % exporting calculated data     
            if ~isempty( obj.outFileName )  
                save( obj.outFileName, 'obj');
            end            
            
        end    
        
        if (idx == obj.MaxIter) && ~obj.ConvergenceTest()
            obj.display( 'EQuUs:SelfConsitent:runSelfConsistentIterations: Iterations failed to converge.') 
        end
        
        
            

        
        
    end    

%% Reset
%> @brief Resets the class members.
    function Reset( obj )

        obj.Initialize();
        
    end
    
    
%% Write
%> @brief Sets the value of an attribute in the class.
%> @param MemberName The name of the attribute to be set.
%> @param input The value to be set.
    function Write(obj, MemberName, input)  
       
        try
          obj.(MemberName) = input;
        catch
               error(['EQuUs:', class(obj), ':Read'], ['No property to write with name: ',  MemberName]);
        end
                
    end
%% Read
%> @brief Query for the value of an attribute in the class.
%> @param MemberName The name of the attribute to be set.
%> @return Returns with the value of the attribute.
    function ret = Read(obj, MemberName)
        
        try
          ret = obj.(MemberName);
        catch
               error(['EQuUs:', class(obj), ':Read'], ['No property to read with name: ',  MemberName]);
        end            
        
    end
%% Clear
%> @brief Clears the value of an attribute in the class.
%> @param MemberName The name of the attribute to be cleared.
    function Clear(obj, MemberName)
          try
          obj.(MemberName) = [];
        catch
               error(['EQuUs:', class(obj), ':Clear'], ['No property to clear with name: ',  MemberName]);
        end  
        
    end        
    
end % public methods


methods ( Access = protected )
    
    
%% Initialize
%> @brief Initializes class attributes.
    function Initialize(obj)
        
        obj.density = [];
        obj.junction_sites = [];  
        obj.filling_factor = [];        
        obj.tolerance = [];
        obj.MaxIter = 1e3;
        
        obj.magnetization = [];
        obj.magnetizaion_direction = [];
        obj.plotIterations = [];
        obj.newIterationWeight = [];
        
        
        iterations_param = 10;
        
        obj.densities_iterations = cell(iterations_param,1);    
        obj.delta_densities_iterations  = nan(iterations_param,1);  
        obj.magnetization_iterations  = nan(iterations_param,1);  
        
        obj.InputParsing( obj.varargin{:} );        
    end            


%% NewIterationInput
%> @brief Creates input for the next self-consistent iteration
        function NewIterationInput(obj)
            
            iterations_param = length(obj.delta_densities_iterations);
            nonzeros = find(obj.delta_densities_iterations); 
            
            if isempty(nonzeros)
                weight_new = 1;
            else            
                weight_new = 0.2;%1/( iterations_param + 1);
                indexes = ~isnan(obj.delta_densities_iterations);
                weights = 1./obj.delta_densities_iterations(indexes);
                if isempty(weights)
                    weights = 1 - weight_new;
                else
                    weights = weights/sum(weights);
                    weights = weights/(sum(weights)+weight_new);
                end
                
            end
            
            
            %> determine new iteration input by weighting the previous results          
            
            obj.density   = weight_new*obj.density;            
            for iidx = 1:length(obj.densities_iterations)
                if isempty( obj.densities_iterations{iidx} )
                    break
                end
                    
                obj.density   =  obj.density + weights(iidx)*obj.densities_iterations{iidx};
            end    
            
            %> renormalizing the density to the initial norm
            filling_factor_tmp = sum(obj.density)/numel(obj.density);
            obj.density   =  obj.density/filling_factor_tmp*obj.filling_factor;
            
            %> calculating the difference compared to the privous density
            if ~isempty( obj.densities_iterations{1} )
                obj.delta_densities_iterations(1) = norm(obj.densities_iterations{1}-obj.density)/norm(obj.densities_iterations{1});
            end
            
            %> calculating the current magnetization
            if obj.Opt.Spin
                scatter_density = obj.density(obj.junction_sites.Scatter.site_indexes);
                obj.magnetization = sum(scatter_density(obj.junction_sites.Scatter.coordinates.spinup) - scatter_density(~obj.junction_sites.Scatter.coordinates.spinup));      
            end
                
            %> shifting iterations factors
            obj.delta_densities_iterations(2:end) = obj.delta_densities_iterations(1:end-1);            
            for iidx = length(obj.densities_iterations):-1:2
                obj.densities_iterations{iidx} = obj.densities_iterations{iidx-1};
            end 
            obj.densities_iterations{1} = obj.density;    
            
            if obj.Opt.Spin
                obj.magnetization_iterations(2:end) = obj.magnetization_iterations(1:end-1);
                obj.magnetization_iterations(1) = obj.magnetization;
            end
            

            
        end


%% ConvergenceTest
%> @brief Tests the convergence of the self-consistent iterations.
%> @return Return with true if convergence was reached, false otherwise
        function ret = ConvergenceTest( obj )
            
            if isnan(obj.delta_densities_iterations(1) )
                ret = false;
                return 
            end
            
            ret1 = obj.delta_densities_iterations(1) < obj.tolerance;   
            
            if ~isempty(obj.Opt.Spin) && obj.Opt.Spin
                if ~isnan( obj.magnetization_iterations(1) ) 
                    ret2 = ( abs((obj.magnetization_iterations(1) - obj.magnetization)/obj.magnetization) < obj.tolerance) ...
                                    || ( abs(obj.magnetization) < obj.tolerance && abs(obj.magnetization_iterations(1)) < obj.tolerance); 
                    ret = ret1 & ret2;
                else
                    ret = false;
                end
                return
            else            
                ret = ret1;
                return;
            end
            
        end

%% input parsing
%> @brief Parses the optional parameters for the class constructor.
%> @param varargin Cell array of optional parameters (https://www.mathworks.com/help/matlab/ref/varargin.html):
%> @param 'resume' Logical value. If true, the iteration continues with the loaded data, or false (default) otherwise.
%> @param 'outFileName' The filename to save the calculated data. If it is empty (default) no data are exported.
%> @param 'silent' Logical value. Set true to suppress the output messages, or false (default) otherwise.
%> @param 'magnetizaion_direction' Set +/- 1 to prescribe the magnetization direction, or set empty (default) to not lock the magnetization direction.
%> @param 'tolerance' Tolerance of the convergence test. (Default is 2*10^-5)
%> @param 'plotIterations' Set true to visuelize the results of the self-consistent iterations.
%> @param 'newIterationWeight' Scalar value. The weith of the newly calculated density in the iterations. (Default value is 0.5)
%> @param 'MaxIter' Maximal number of iterations. (Default value is 500)
    function InputParsing( obj, varargin )
        
        p = inputParser;
        p.addParameter('resume', 0);
        p.addParameter('outFileName', []);
        p.addParameter('silent', obj.Opt.Silent);
        p.addParameter('magnetizaion_direction', []);
        p.addParameter('tolerance', 2e-5);
        p.addParameter('plotIterations', false);
        p.addParameter('newIterationWeight', 0.5);
        p.addParameter('MaxIter', 500);
        
        p.parse(varargin{:});

       
        obj.resume        = p.Results.resume;
        obj.outFileName       = p.Results.outFileName;
        obj.silent  = p.Results.silent;
        obj.magnetizaion_direction = p.Results.magnetizaion_direction;
        obj.tolerance = p.Results.tolerance;
        obj.plotIterations = p.Results.plotIterations;
        obj.newIterationWeight = p.Results.newIterationWeight;
        obj.MaxIter = p.Results.MaxIter;
        
        
    end



end % methods private


end