Homec4science
Diffusion lammps (lammm-master) Browse

lammps/tools/matlab/readEAM.m
0f1efa0c79d6lammm-master

Last Change
Raw File

readEAM.m
View Options

function varargout = readEAM(varargin)
% Function to read EAM potential files
% Input
% 1: EAM potential file name
% 2: file type --> 'FUNCFL' for single element file
% --> 'SETFL' for multiple element file
% Output : Structure with members
% .embed : Embedding function
% .pair : Pair Potential
% .elecden : Electron Density
% .nrho : Number of points for electron density
% .drho : increment of r for electron density
% .nr : Number of points for pair potential
% .dr : increment of r for pair potential
% .rcut : cut-off distance
% All output is in exactly the same units as in the EAM file. For
% multielement SETFL files the members are multidimensional arrays with
% each element data stored in (:,:,ELEM)
%
% Example
% eam = readEAM('cuu3.eam','funcfl');
%
% Author : Arun K. Subramaniyan
% sarunkarthi@gmail.com
% http://web.ics.purdue.edu/~asubrama/pages/Research_Main.htm
% School of Aeronautics and Astronautics
% Purdue University, West Lafayette, IN - 47907, USA.
if length(varargin) < 2
error('Too few input parameters : Need Filename and filetype');
elseif length(varargin) > 2
error('Too many input parameters');
end
filename = varargin{1};
type = varargin{2};
try
fid = fopen(filename,'r');
catch
error('EAM file not found!');
end
if strcmpi(type,'FUNCFL')
t = 1;
elseif strcmpi(type,'SETFL')
t = 2;
else
error(['Unknown file type : "' type '"']);
end
switch t
case 1 %FUNCFL
% Read line 1 : comment line
a = fgetl(fid);
% read Line 2
rem = fgetl(fid); % ielem amass blat lat
[ielem,rem] = strtok(rem);
[amass,rem] = strtok(rem);
[blat,rem] = strtok(rem);
[lat,rem] = strtok(rem);
ielem = str2num(ielem); % Atomic Number
amass = str2num(amass); % Atomic Mass
blat = str2num(blat); % Lattice Constant
% Read Line 3
a = str2num(fgetl(fid));
nrho = a(1);
drho = a(2);
nr = a(3);
dr = a(4);
rcut = a(5);
% Reading embedding function
for i = 1 : 1 : nrho/5
embed(i,:) = str2num(fgetl(fid));
end
% Reading pair potential
for i = 1 : 1 : nr/5
pair(i,:) = str2num(fgetl(fid));
end
% Reading electron density function
for i = 1 : 1 : nr/5
elecden(i,:) = str2num(fgetl(fid));
end
% Output
out.embed = embed;
out.pair = pair;
out.elecden = elecden;
out.nrho = nrho;
out.drho = drho;
out.nr = nr;
out.dr = dr;
out.rcut = rcut;
varargout{1} = out;
% --------------------------------------------------------------------
case 2 % SETFL
% Read lines 1 - 3 : comment lines
a = fgetl(fid);
a = fgetl(fid);
a = fgetl(fid);
% Atom types
ntypes = str2num(fgetl(fid));
% Read Global information
a = str2num(fgetl(fid));
nrho = a(1);
drho = a(2);
nr = a(3);
dr = a(4);
rcut = a(5);
% Read element specific Data
% Embedding function and Electron Density
for elem = 1 : 1 : ntypes
rem = fgetl(fid); % ielem amass blat lat
[ielem1,rem] = strtok(rem);
[amass1,rem] = strtok(rem);
[blat1,rem] = strtok(rem);
[lat1,rem] = strtok(rem);
ielem(elem) = str2num(ielem1); % Atomic Number
amass(elem) = str2num(amass1); % Atomic Mass
blat(elem) = str2num(blat1); % Lattice Constant
lat(elem,:) = lat1; % Lattice type
% Reading embedding function
for i = 1 : 1 : nrho/5
embed(i,:,elem) = str2num(fgetl(fid));
end
% Reading electron density function
for i = 1 : 1 : nr/5
elecden(i,:,elem) = str2num(fgetl(fid));
end
end
% Pair Potentials
n_pair = ntypes * (ntypes + 1) / 2;
for np = 1 : 1 : n_pair
for i = 1 : 1 : nr/5
pair(i,:,np) = str2num(fgetl(fid));
end
end
% Output
out.embed = embed;
out.elecden = elecden;
out.pair = pair;
out.nrho = nrho;
out.drho = drho;
out.nr = nr;
out.dr = dr;
out.rcut = rcut;
varargout{1} = out;
end
c4science · Help