cdse-analysis-scripts/process_log.m

#!/usr/bin/env octave

% Copyright (C) 2019 Kei Kebreau
%
% 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/>.

% Let user select input file.
[input_file, file_path] = uigetfile({'*.log*', 'Supported file formats'});
% Load the input file.
raw_data = fileread([file_path input_file]);
[~, file_name, ~] = fileparts(input_file);
% Let user select output directory.
output_dir = uigetdir('.', 'Select output directory for graph and xyz files');

% Get the number of atoms.
num_atoms = regexp(raw_data, 'NAtoms=\s+(\d+)', 'tokens', 'once');
% Extract the string from the cell array and convert it to a number.
num_atoms = str2num(num_atoms{1});

% There *has* to be a better way to do this, but for now it works... :-/
intermediate_cluster_regexp = ['Standard orientation:' repmat('\D+\d+\s+(\d+)\s+\d+\s+([+-]?[0-9]*\.?[0-9]+)\s+([+-]?[0-9]*\.?[0-9]+)\s+([+-]?[0-9]*\.?[0-9]+)', [1, num_atoms])];
intermediate_clusters = regexp(raw_data, intermediate_cluster_regexp, 'tokens');
intermediate_clusters = cellfun(@str2num, [intermediate_clusters{:}]);
steps = length(intermediate_clusters) / num_atoms / 4;

% Also determine how many leading zeroes are necessary for correctly listing
% the intermediate .xyz files.
leading_zeroes = floor(log10(steps));

previous_energy = 0;
step_count = 0;
xc_functional = regexp(raw_data, 'SCF Done:\s+E\((\w+)\)', 'tokens', 'once');
xc_functional = xc_functional{1};
% Extract the days, hours, minutes and seconds.
cpu_time = regexp(raw_data,
                  'Job cpu time:\D+(\d+)\D+(\d+)\D+(\d+)\D+([0-9]*\.?[0-9]+)',
                  'tokens', 'once');
% Convert to 'days:hours:minutes:seconds' format.
cpu_time = [cpu_time{1} ':' cpu_time{2} ':' cpu_time{3} ':' cpu_time{4}];
energy = regexp(raw_data, 'SCF Done:\s+E\(\w+\) =  ([+-]?[0-9]*\.?[0-9]+)', 'tokens');
energy = cellfun(@str2num, [energy{:}]);
pct_change_in_energy = -diff(energy) ./ energy(1:end-1) * 100;
cycles = regexp(raw_data, 'SCF Done:\s+\S+\s+\S+\s+\S+\s+\S+\s+\S+\s+(\d+)', 'tokens');
cycles = cellfun(@str2num, [cycles{:}]);

% Remember: intermediate_clusters has all atomic numbers and 3D locations in a flat matrix!
% Parse wisely.
xyz_format_str = ['%s%s%s_Step_%0' num2str(leading_zeroes) 'd.xyz'];
for step = 1:steps
  xyz_name = sprintf(xyz_format_str, output_dir, filesep, file_name, step);
  xyz_id = fopen(xyz_name, 'w');
  fdisp(xyz_id, sprintf('%d\n', num_atoms));
  for atom = 1:num_atoms
    fdisp(xyz_id, sprintf('%d %f %f %f',
                          intermediate_clusters((step - 1) * num_atoms * 4
                                                + (atom - 1) * 4
                                                + 1),
                          intermediate_clusters((step - 1) * num_atoms * 4
                                                + (atom - 1) * 4
                                                + 2),
                          intermediate_clusters((step - 1) * num_atoms * 4
                                                + (atom - 1) * 4
                                                + 3),
                          intermediate_clusters((step - 1) * num_atoms * 4
                                                + (atom - 1) * 4
                                                + 4)));
  endfor
  fclose(xyz_id);
endfor

figure(1, 'position', get(0, 'screensize'), 'visible', 'off');
subplot(2, 1, 1);
plot(cycles, energy, 'o');
xlabel('Cycles');
ylabel('Energy (A.U.)');
title(['Cycles at Different Energy Levels (Functional: ' xc_functional ')']);
set(gca, 'fontsize', 16);
set(gca, 'xdir', 'reverse');
subplot(2, 1, 2);
plot(cycles(2:end), pct_change_in_energy, 'o');
xlabel('Cycles');
ylabel('% Change in Energy (A.U.)');
set(gca, 'fontsize', 16);
set(gca, 'xdir', 'reverse');
print([output_dir filesep file_name '_energy_graph.png']);
hgsave([output_dir filesep file_name '_energy_graph.ofig']);
set(gcf, 'visible', 'on');