cdse-analysis-scripts/generate_cdse.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/>.

% Return a matrix representation of a tetrahedral CdSe cluster with a given
% number of layers of cadmium, and optionally save the cluster as an XYZ file
% specified by the given file name.
function cdse_cluster = generate_cdse (layer_count, file_name)
  try
    pkg load geometry;
  catch
    error('Could not load Octave Forge package "geometry"');
  end
  % The number of Cd atoms in a cluster of 'n' layers is the cumulative sum of
  % the first 'n' triangular number for all 'n' > 0 while the number of Se atoms
  % in a cluster of 'n' layers is the cumulative sum of the first 'n-1'
  % triangular numbers for all 'n' > 1.
  if layer_count < 2
    error('A tetrahedral CdSe cluster requires at least 2 layers.');
  end
  se_count = sum(arrayfun(@triangular_number, 1:layer_count-1));
  cd_count = se_count + triangular_number(layer_count);
  % Create a tetrahedral grid.
  [vertices, edges, faces] = createTetrahedron;
  cd_scaling_factor = 3.0385 * (layer_count - 1);
  se_scaling_factor = 3.0385 * (layer_count - 2);
  cd_cluster = tetrahedron_grid(layer_count - 1, vertices', cd_count)';
  se_cluster = tetrahedron_grid(layer_count - 2, vertices', se_count)';
  cd_cluster = cd_cluster * cd_scaling_factor;
  se_cluster = se_cluster * se_scaling_factor + [1.5193 1.5193 1.5193];
  cdse_cluster = [cd_cluster; se_cluster];

  % Optional file output.
  if nargin < 2
    return;
  end
  xyz_id = fopen(file_name, 'w');
  fdisp(xyz_id, sprintf('%d\n', cd_count + se_count));
  for ii = 1:size(cd_cluster)(1)
    fdisp(xyz_id, sprintf('Cd %f %f %f', cd_cluster(ii,:)));
  end
  for ii = 1:size(se_cluster)(1)
    fdisp(xyz_id, sprintf('Se %f %f %f', se_cluster(ii,:)));
  end
  fclose(xyz_id);
end

% Return the nth triangular number.
function result = triangular_number (n)
  result = n * (n + 1) / 2;
end