function [ unit_list,all_e2_sum_matrix,all_e2_max_matrix,all_e2_max_info,h_e2_sum_matrix,h_e2_max_matrix,h_e2_max_info ] = NBO_E2_Vals_File( filename )
%Function extracts E2 energy values (in kJ mol) for NBO.log files
%and puts them in a more useful format (also outputs sums of E2 vals).

%********************************Summary********************************%
%   Extracts useful info from the 2nd order pertubation analysis        %
%	section of an NBO .log file. Gives outputs which are the sum        %
%   of all D/A interactions between units, as well as the strongest     %
%   D/A interaction and a list of all D/A interactions. Furthermore     %
%   outputs the same info but only taking D/A interactions where        %
%   MOs with Hydrogen contribs act as the acceptor(i.e as they do in    %
%   hydrogen bonds).                                                    %
%*********************************Input*********************************%
%   1) filename = strnig containing full path to the NBO file of        %
%                 interest.                                             %
%********************************Output*********************************%
%   1) unit_list = nx1 (column) cell array containing names of all      %
%                   "units" (should be molecules) in the .log file      %
%   2) all_e2_sum_matrix = NxN matrix (n=number of units in file).      %
%                          Organised as follows(LHS=donor):             %
%       [unit 1--> unit 1][unit 1--->unit 2]                            %
%       [unit 2--> unit 1][unit 2--->unit 2]                            %
%                          Each element contains the sum of ALL the     %
%                          relevant D-->A e2 energies in kJ mol^-1      %
%                          units of kJ mol ^{-1}                        %
%   3) all_e2_max_matrix = same as all_e2_sum matrix, except matrix     %   
%                          elements correspond to the highest energy    %
%                          interaction rather than the sum              %
%                          units of kJ mol ^{-1}                        %
%   4) all_e2_max_info = cell array (NxN,N=number units) in same        % 
%                        format as all_e2_max_matrix, but each element  %
%                        contains the entirerity of the output line     %
%                        units of kcal mol ^{-1}                        %
%   5) h_e2_sum_matrix = same as all_e2 except with just interactions   %
%                        where H atoms act as acceptor being considered %
%   6) h_e2_max_matrix = same as all_e2 version, except only            %
%                        interactions with H atoms act as acceptor      %
%                        units of kJ mol ^{-1}                          %
%   7) h_e2_max_info = same as all_e2 version, except only              %
%                      interactions with H atoms act as acceptor        %
%                      units of kcal mol^{-1}                           %
%********************************Issues*********************************%
%********************************History********************************%
%***********************************************************************%






    %********************************Summary********************************%
    %   Extracts useful info from the 2nd order pertubation analysis        %
    %	section of an NBO .log file (also gets definitions of "units"  	    %
    %	used in the pert. section.                                          %
    %*********************************Input*********************************%
    %	1) filename = string containing the full path to the NBO .log file  %
    %*********************************Output********************************%
    %   1) unit_list = 1 x n (column) cell array containing mol. formula    %
    %                  of units 1-->n                                       %
    %   2) full_line = m x (n^2) cell array (n=number of units in file,     %
    %                  m is number of donor/acceptor interactions between   %
    %                  any two units, and hence unknown in advance.         %
    %		       Each entry is a full line from the E(2) pert.            %
    %		       section of the NBO file. column 1 is 1-->1 CT,           %
    %		       2nd is 1-->2 CT, nth is 1-->n CT, n+1 col is 2-->1   	%
    %		       etc. etc.                                                %
    %   3) e2_vals = m x (n^2) matrix, same format as full_line but         %
    %		     with only the E2 values (in kJ mol)                        %
    %   4) h_full_line = same as full_line except only interactions with    %
    %                    hydrogen acting as a donor are printed             %
    %   5) h_e2_vals = same as e2 vals, but with only e2 values for         %
    %                  D/A interactions with H as an acceptor               %
    %***********************************************************************% 
    %*********************************Issues********************************%
    %   1) Suspect that it may fail for numb_units>10 (due to the           %
    %      search_string_matrix not working. Suspect it needs to be         %
    %      changed to search for a reg. expr. to solve                      %
    %***********************************************************************% 

    function [unit_list,full_line,e2_vals,h_full_line,h_e2_vals] = import_file_NBO(filename)

        %****SECTION 1****%
        %Put file in 1xn cell array(column), each line=1 row
        fobject1=fopen(filename);
        full_file = fileread(filename);
        full_file = strsplit(full_file,'\n')';
        fclose(fobject1);
        
        %****SECTION 2****%
        %Get line numbers for all lines which define chem. formula of a "unit"
        lnumb_string = ''; %comma delim string
        search_term = 'Molecular unit' ;
        for i=1:size(full_file,1)
            if ~isempty(strfind(full_file{i,1},search_term));
                if strcmp(lnumb_string,'')
                    lnumb_string = num2str(i);
                else
                    lnumb_string = strcat(lnumb_string,',',num2str(i));
                end
            end
        end
        
        %****SECTION 3****%
        %Extract chemical formula for each unit
        lnumb_array = strsplit(lnumb_string,',')';	%1xn (column) cell array w relevant line numbers
        unit_list = cell(size(lnumb_array,1),1) ; %OUTPUT ARG. see funct. description
        for i=1:size(lnumb_array,1)
            curr_line = str2num(lnumb_array{i,1});
            unit_list{i,1} =regexp(full_file{curr_line,1},'\([a-za-zA-Z_0-9]*)','match');
        end
        
        
        %****SECTION 4****%
        %Cut file down to the section we need.(store in cotnracted_file cell array)
        for i = 1 :size(full_file,1)
            if ~isempty(strfind(upper(full_file{i,1}),'SECOND ORDER PERTURBATION'));
                start_line = i+6;
                for j=i+1:size(full_file,1)
                    if  ~isempty(strfind(upper(full_file{j,1}),'NATURAL BOND ORBITALS'));
                        end_line = j-1;
                        break;
                    end
                end
                break;
            end
        end
        contracted_file = full_file(start_line:end_line,1);
        
        %****SECTION 5****%
        %Get full_line output variable
        %NOTE: search_string and line number marices have format:
        %[unit 1 -->unit 1]
        %[unit 1-->unit 2]
        %[unit 1-->unit n]
        %[unit 2--->unit 1]
        %[unit 2-->unit 2] etc.
        numb_units = size(unit_list,1);
        search_string_matrix = cell(numb_units*numb_units,1); %nx1 (column) cell array containg a string for each type of unit donor/acceptor interaction
        line_number_matrix = NaN(numb_units*numb_units,2) %nx2 (column) matrix, 1st col = start line (in contreaced_file) for the info for unit in search_string, 2nd=end line
        line_number_matrix(end,2) = size(contracted_file,1);
        %Generate search_string
        counter=1;
        for donor=1:numb_units
            for acceptor=1:numb_units
                if donor==acceptor
                    search_string_matrix{counter,1} = ['within unit  ' num2str(acceptor)];
                    counter = counter+1;
                else
                    search_string_matrix{counter,1} = ['from unit  ' num2str(donor) ' to unit  ' num2str(acceptor)];
                    counter=counter+1;
                end
            end
        end
        
        %Find start/end points of each donor/acceptor relationship
        %-------------->NOTE:a val of 0 in line_number_matrix(x,2) means no CT for that d/A relationship<---------------%
        counter=1;
        for i=1:size(search_string_matrix,1)
            instances_found = 0;		%To make sure the strings only found once (Error checking)
            for j=1:size(contracted_file,1)
                if strcmp([' ' search_string_matrix{i,1}],lower(contracted_file{j,1}))
                    if instances_found>0
                        error(['Multiple donor/acceptor relationships found for string '  search_string_matrix{i,1}]);
                    else
                        
                        instances_found=instances_found+1;
                        line_number_matrix(i,1) = j+1;	%line where this d/a relatinoship info starts (in contr. file)
                        %Next for loop is to find the end positions
                        for k=j+1:size(contracted_file,1);
                            if isempty(regexp(contracted_file{k,1},'\d\.'))
                                line_number_matrix(i,2) = k-1;	% final line containing values
                                if strfind(contracted_file{k,1},'None') %Test for case of no CT above threshold
                                    line_number_matrix(i,2)=0;
                                end
                                break;
                            end
                        end
                    end
                end
            end
        end
        
        %Put relevant lines in full_lines variable
        full_line = cell(max((line_number_matrix(:,2) - line_number_matrix(:,1))),numb_units*numb_units);
        for i=1:size(line_number_matrix,1)
            temp_array = contracted_file(line_number_matrix(i,1):line_number_matrix(i,2),1);
            full_line(1:size(temp_array,1),i) = temp_array;
        end
        
        
        %****SECTION 6****%
        %Get e2_vals/h_full_line/h_e2_vals
        
        e2_vals = NaN(size(full_line,1),size(full_line,2));
        h_e2_vals = NaN(size(full_line,1),size(full_line,2));
        h_full_line = cell(size(full_line,1),size(full_line,2));
        
        for column = 1:size(full_line,2)
            h_counter = 1;	%counts the row we're on in the h_{e2_vals/full_line}
            for i=1:size(full_line,1)
                curr_line = full_line{i,column};
                if ~isempty(curr_line)
                    donor = curr_line(1:22);	%str containing info on donor MO
                    acceptor = curr_line(29:51);			%str contaninig info on acceptr MO
                    numbers = curr_line(55:end);
                    numbers_array = strsplit(strtrim(numbers),' '); %1st col=E(2),2nd col = ej-ei, 3rd col = fij
                    e2_vals(i,column) = str2num(numbers_array{1,1});
                    %Now extracting e2 val/line if its got a hydrogen acceptor in it
                    if ~isempty(strfind(acceptor,'H'))
                        h_e2_vals(h_counter,column) = str2num(numbers_array{1,1});
                        h_full_line{h_counter,column} = curr_line;
                        h_counter=h_counter+1;
                    end
                end
            end
        end
        %converting to kj mol (from kcal)
        e2_vals = e2_vals *4.1840;
        h_e2_vals = h_e2_vals *4.1840;
        
        %Deletes the rows that are totally empty in h_{e2_vls/full_line}
        max_rows = 0;	%max number of non-empty rows in the h_* matrix/array
        for j=1:size(h_full_line,2);
            for i=1:size(h_full_line,1);
                if isempty(h_full_line{i,j});
                    if (i-1) > max_rows
                        max_rows = i-1;
                    end
                    break	;
                end
            end
        end
        h_full_line = h_full_line (1:max_rows,:);
        h_e2_vals = h_e2_vals (1:max_rows,:);
        
    end


    %********************************Summary********************************%
    %   Extracts the information on the strongest donor/acceptor            %
    %   interaction between all units                                       %
    %*********************************Input*********************************%
    %   1)e2_vals = mxn matrix containing e2 values (and NaN if no val      %
    %               available). same as output from "import_file_NBO".      %
    %               1 col per D/A unit interaction.column 1 is 1-->1 CT,    %
    %		       2nd is 1-->2 CT, nth is 1-->n CT, n+1 col is 2-->1   	%
    %		       etc. etc.                                                %
    %   2)e2_max_info = mxn cell array, each entry containing the full      %
    %                   line in the .log file                               %
    %*********************************Output********************************%
    %   1)e2_max_matrix = matrix containing the maximum e2 values between   %
    %                     units in the format below:                        %
    %                     [unit 1-->unit 1] [unit 1-->unit 2]               %
    %                     [unit 2-->unit 1] [unit 2-->unit 2]               %
    %                     E(2) vals here are kJ mol^{-1}
    %   2)e2_max_info = cell array in same format as e2_max_matrix,         %
    %                   except each entry contains the full line from       %
    %                   the input file.NOTE E(2) vals in this are in kcal   %
    %                   per mol
    %***********************************************************************% 
    function [e2_max_matrix,e2_max_info] = gen_max_e2_info(e2_vals,e2_info)
        
        %need to get number of frags 1st:
        numb_frags = sqrt(size(e2_info,2));
        if rem(numb_frags,1) ~= 0 ;
            error('Invalid number of fragments');
        end
        
        %Extracting max E2 value for each donor.acceptor relationship
        e2_max_matrix = NaN(numb_frags,numb_frags);
        e2_max_info = cell(numb_frags,numb_frags);
        counter=1;	%tracks column of e2_vals/info we're on
        %row and col loop vars are the row/col of the output vars
        for row = 1:numb_frags
            for col=1:numb_frags
                [max_e2_val,loc] = max(e2_vals(:,counter))	%loc is the row that contains the e2 val
                e2_max_matrix(row,col) = e2_vals(loc,counter);
                e2_max_info(row,col) =  e2_info(loc,counter);
                counter = counter+1;
            end
        end
    end



    %********************************Summary********************************%
    %   Extracts the information on the sum of E(2) values for each         %
    %   different combination of units                                      %
    %*********************************Input*********************************%
    %   1)e2_vals = mxn matrix containing e2 values (and NaN if no val      %
    %               available). same as output from "import_file_NBO".      %
    %               1 col per D/A unit interaction.column 1 is 1-->1 CT,    %
    %		       2nd is 1-->2 CT, nth is 1-->n CT, n+1 col is 2-->1   	%
    %		       etc. etc.                                                %
    %*********************************Output********************************%
    %   1)e2_sum_matrix = matrix containing the summed e2 values between    %
    %                     units in the format below:                        %
    %                     [unit 1-->unit 1] [unit 1-->unit 2]               %
    %                     [unit 2-->unit 1] [unit 2-->unit 2]               %
    %                     E(2) vals here are kJ mol^{-1}                    %
    %***********************************************************************% 
    function [e2_sum_matrix] = gen_sum_e2_matrix(e2_vals)
        %need to get number of frags 1st:
        numb_frags = sqrt(size(e2_vals,2));
        if rem(numb_frags,1) ~= 0 ;
            error('Invalid number of fragments');
        end
        
        
        %Extracting max E2 value for each donor.acceptor relationship
        e2_sum_matrix = NaN(numb_frags,numb_frags);
        
        counter=1;	%tracks column of e2_vals/info we're on
        %row and col loop vars are the row/col of the output vars
        for row = 1:numb_frags
            for col=1:numb_frags
                e2_sum_matrix(row,col) = nansum(e2_vals(:,counter));
                counter = counter+1;
            end
        end
    end

    %obv see top of file for description
    function [] = Main()
           [unit_list,full_line,e2_vals,h_full_line,h_e2_vals] = import_file_NBO(filename);
           [all_e2_sum_matrix] = gen_sum_e2_matrix(e2_vals);
           [all_e2_max_matrix,all_e2_max_info] = gen_max_e2_info(e2_vals,full_line);
           [h_e2_sum_matrix] = gen_sum_e2_matrix(h_e2_vals);
           [h_e2_max_matrix  h_e2_max_info] = gen_max_e2_info(h_e2_vals,h_full_line);
    end


tic
Main
toc


end