diff --git a/tools/matlab/README.pdf b/tools/matlab/README.pdf new file mode 100644 index 0000000000..e10c05a5a0 Binary files /dev/null and b/tools/matlab/README.pdf differ diff --git a/tools/matlab/lmp2cfg.m b/tools/matlab/lmp2cfg.m new file mode 100644 index 0000000000..743fd90c22 --- /dev/null +++ b/tools/matlab/lmp2cfg.m @@ -0,0 +1,252 @@ +function lmp2cfg(varargin) +% Converts LAMMPS dump file to Extended CFG Format (No velocity) to be used +% with AtomEye (http://164.107.79.177/Archive/Graphics/A/) +% Input : +% Necessary (in order) +% timestep,Natoms,x_bound,y_bound,z_bound,H,atom_data,mass, +% cfg file name, dumpfile name +% Optional +% 'vel' --> 'yes' (if velocity data needs to be written) +% --> Default is 'no' +% 'atomtype' --> {Atom symbol} +% 'autonumber' --> 'yes' (default) numbers cfgfiles consecutively based on timestep +% 'autolog' --> 'on' (default) writes details of conversion to a log file +% 'rotation' --> To specify rotation (Transform in cfg file) say 'rotation' followed by +% number of rotations,the axes of rotation ('X' or 'Y' or 'Z')and the +% angle in degrees. CAUTION : Do remember that rotations are +% non-commutative and hence the order of rotations must be the same as +% intended +% 'aux' --> Auxiliary data must be a two column array with Quantity and Type in each +% column +% +% +% THE DATA MUST BE SCALED (0 to 1)IN ORDER FOR ATOMEYE TO READ PROPERLY +% Real coordinates x = s * H, x, s are 1x3 row vectors +% +% Example +% lmp2cfg_all(data.timestep,data.Natoms,data.xbound,data.ybound,... +% data.zbound,H,data.atom_data,mass,cfgfile,dumpfile,... +% 'vel','no','aux',{'sxx' 'Lammps O/P' +% 'syy' 'Lammps O/P' +% 'szz' 'Lammps O/P' +% 'sxy' 'Lammps O/P' +% 'sxz' 'Lammps O/P' +% 'syz' 'Lammps O/P'},... +% 'atomtype',{'Ni' +% 'Al'},... +% 'autonumber','on',... +% 'autolog','on',... +% 'rotation',2,'X',45,'Y',30); +% +% See also readdump_all, readdump_one, scandump +% +% 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. + +%------------ Defaults +vel = 'no'; +auxstatus = 0; +atom_typestatus = 0; +rotstatus = 0; +autonumber_status = 1; % default is ON +autolog_status = 1; +%---------------------------------- + +%-----Required Input -------------- +timestep = varargin{1} +Natoms = varargin{2}; +x_bound = varargin{3}; +y_bound = varargin{4}; +z_bound = varargin{5}; +H = varargin{6}; +atom_data = varargin{7}; +mass = varargin{8}; %arrange it in order of atom type , ie 1 - mass1 etc +filename = varargin{9}; %Just give one file name + +if length(varargin) < 10 + dumpfilename = ['dump.' filename]; +else + dumpfilename = varargin{10}; +end + + + +if length(varargin) > 10 + i=11; + while i \t\t %s \n',timestep(i),name{i}); + end + fclose(flog); +end + +%---------- Function to calculate beta +function b = beta(angle,axes) + switch axes + case 'X' % X axes + b = [1 0 0 + 0 cos(angle) sin(angle) + 0 -sin(angle) cos(angle)]; + case 'Y' % Y axes + b = [cos(angle) 0 sin(angle) + 0 1 0 + -sin(angle) 0 cos(angle)]; + case 'Z' % Z axes + b = [cos(angle) sin(angle) 0 + -sin(angle) cos(angle) 0 + 0 0 1]; + end +end +% -------------------------------------------------- + +end % For main function + + diff --git a/tools/matlab/readEAM.m b/tools/matlab/readEAM.m new file mode 100644 index 0000000000..92c39a9629 --- /dev/null +++ b/tools/matlab/readEAM.m @@ -0,0 +1,166 @@ +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 + (factorial(ntypes)/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 + diff --git a/tools/matlab/readdump_all.m b/tools/matlab/readdump_all.m new file mode 100644 index 0000000000..4f8270b4f0 --- /dev/null +++ b/tools/matlab/readdump_all.m @@ -0,0 +1 @@ +function [varargout] = readdump_all(varargin) % Reads all timesteps from a LAMMPS dump file. % Input is dump file name with path % Output is in the form of a structure with following variables % .timestep --> Vector containing all time steps % .Natoms --> Vector containing number of atoms at each time step % .x_bound --> [t,2] array with xlo,xhi at each time step % .y_bound --> [t,2] array with ylo,yhi at each time step % .z_bound --> [t,2] array with zlo,zhi at each time step % .atom_data --> 3 dimensional array with data at each time step stored % as atomdata(:,:,t) % Example % data = readdump_all('dump.LAMMPS'); % % See also readdump_one, scandump % % 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. try dump = fopen(varargin{1},'r'); catch error('Dumpfile not found!'); end i=1; while feof(dump) == 0 id = fgetl(dump); switch id case 'ITEM: TIMESTEP' timestep(i) = str2num(fgetl(dump)); case 'ITEM: NUMBER OF ATOMS' Natoms(i) = str2num(fgetl(dump)); case 'ITEM: BOX BOUNDS' x_bound(i,:) = str2num(fgetl(dump)); y_bound(i,:) = str2num(fgetl(dump)); z_bound(i,:) = str2num(fgetl(dump)); case 'ITEM: ATOMS' for j = 1 : 1: Natoms atom_data(j,:,i) = str2num(fgetl(dump)); end i=i+1; end end %----------Outputs------------- %OUTPUTS IN SAME VARIABLE STRUCTURE varargout{1}.timestep = timestep; varargout{1}.Natoms = Natoms; varargout{1}.x_bound = x_bound; varargout{1}.y_bound = y_bound; varargout{1}.z_bound = z_bound; varargout{1}.atom_data = atom_data; %------------------------------ fclose(dump); \ No newline at end of file diff --git a/tools/matlab/readdump_one.m b/tools/matlab/readdump_one.m new file mode 100644 index 0000000000..f5b1a13354 --- /dev/null +++ b/tools/matlab/readdump_one.m @@ -0,0 +1,120 @@ +function [varargout] = readdump_one(varargin) +% Read LAMMPS dump file one timestep at a time +% Input +% Dump file name with path +% Starting file pointer position in dump file +% Number of columns in the dump file +% Output is in the form of a structure with following variables +% .timestep --> Vector containing all time steps +% .Natoms --> Vector containing number of atoms at each time step +% .x_bound --> [t,2] array with xlo,xhi at each time step +% .y_bound --> [t,2] array with ylo,yhi at each time step +% .z_bound --> [t,2] array with zlo,zhi at each time step +% .atom_data --> 2 dimensional array with data +% .position --> file pointer for reading next time +% +% Example +% data = readdump_one('dump.LAMMPS',0,5); +% Reads the first timestep in the file dump.LAMMPS +% Specify position = -1 if only the last dump step in the +% file is needed +% +% See also readdump, scandump +% +% 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. + +try + dump = fopen(varargin{1},'r'); +catch + error('Dumpfile not found!'); +end +position = varargin{2}; % from beg of file +ncol = varargin{3}; %number of columns + +i=1; +t=0; +done = 0; +last_status = 0; +if position ~= -1 + fseek(dump,position,'bof'); +else + last_status = 1; +end +while done == 0 & last_status == 0 + id = fgetl(dump); + switch id + case 'ITEM: TIMESTEP' + if t == 0 + timestep(i) = str2num(fgetl(dump)); + t=1; + end + case 'ITEM: NUMBER OF ATOMS' + Natoms = str2num(fgetl(dump)); + case 'ITEM: BOX BOUNDS' + x_bound(1,:) = str2num(fgetl(dump)); + y_bound(1,:) = str2num(fgetl(dump)); + z_bound(1,:) = str2num(fgetl(dump)); + case 'ITEM: ATOMS' + atom_data = zeros(Natoms,ncol);%Allocate memory for atom data + for j = 1 : 1: Natoms + atom_data(j,:) = str2num(fgetl(dump)); + end + done = 1; + p = ftell(dump); + end +end + +% Getting only the last step +if last_status == 1 + % First get the position of the beginning of the last step in the + % dumpfile + while ~feof(dump) + temp = fgetl(dump); + if length(temp) == 14 + if strcmpi(temp,'ITEM: TIMESTEP') + p = ftell(dump); % starting position of line next to the header + end + end + end + fclose(dump); + dump = fopen(varargin{1},'r'); + fseek(dump,p,'bof'); + % getting Timestep + timestep = str2num(fgetl(dump)); + + while ~feof(dump) + id = fgetl(dump); + switch id + case 'ITEM: NUMBER OF ATOMS' + Natoms = str2num(fgetl(dump)); + case 'ITEM: BOX BOUNDS' + x_bound(1,:) = str2num(fgetl(dump)); + y_bound(1,:) = str2num(fgetl(dump)); + z_bound(1,:) = str2num(fgetl(dump)); + case 'ITEM: ATOMS' + atom_data = zeros(Natoms,ncol);%Allocate memory for atom data + for j = 1 : 1: Natoms + atom_data(j,:) = str2num(fgetl(dump)); + end + end + end +end + +%----------Outputs------------- + +%OUTPUTS IN SAME VARIABLE STRUCTURE +varargout{1}.timestep = timestep; +varargout{1}.Natoms = Natoms; +varargout{1}.x_bound = x_bound; +varargout{1}.y_bound = y_bound; +varargout{1}.z_bound = z_bound; +varargout{1}.atom_data = atom_data; +varargout{1}.position = p; %gives postion of ITEM: TIMESTEP line +%------------------------------ + +fclose(dump); + diff --git a/tools/matlab/readlog.m b/tools/matlab/readlog.m new file mode 100644 index 0000000000..3392d3b267 --- /dev/null +++ b/tools/matlab/readlog.m @@ -0,0 +1,92 @@ +function [varargout] = readlog(varargin) +% Read LAMMPS log files +% input is log file name with path +% output is a structure --> out +% out.Chead --> has heading of columns in each run +% out.data --> has the data during each run in the form of string +% Type str2num(out.data{i}) to get the numeric array +% +% Example +% logdata = readlog('log.LAMMPS'); +% +% Author : sarunkarthi@gmail.com +% http://web.ics.purdue.edu/~asubrama/pages/Research_Main.htm +% Arun K. Subramaniyan +% School of Aeronautics and Astronautics +% Purdue University, West Lafayette, IN - 47907, USA. + +logfile = varargin{1}; +try + fid = fopen(logfile,'r'); +catch + error('Log file not found!'); +end + +loop = 1; +while feof(fid) == 0 + %----------- To get first line of thermo output -------- + while feof(fid) == 0 + a = fgetl(fid); + if length(a) > 4 && strcmp(a(1:4),'Step') + findstep = 1; + break; + end + end + %------------------------------------------------------- + + %---------Seperate column headings---------------------- + j=1; + k=1; + Ch=''; + i=1; + while i < length(a) && findstep == 1 + for i = k : 1 : length(a) + if strcmp(a(i),' ') + Chead{loop,j} = Ch; + clear Ch; + Ch = ''; + k=i+1; + j=j+1; + break; + else + Ch = [Ch a(i)]; + end + end + end + if findstep == 1 + Chead{loop,j} = Ch; % to get the last column heading + end + %------------------------------------------------------- + + %----------------------Get Data------------------------- + id = 1; %row id... + while feof(fid) == 0 + a = fgetl(fid); + if strcmp(a(1:4),'Loop') + loop = loop + 1; + break; + else + logdata(id,:) = str2num(a); + id = id+1; + end + end + %-------------------------------------------------------- + if feof(fid) == 0 + b = num2str(logdata); + data{loop-1} = b; + clear b; + clear logdata; + findstep = 0; + end +end +fclose(fid); + +%--------OUTPUT------------------------------------------- +out.Chead = Chead; +out.data = data; + +varargout{1} = out; + + + + \ No newline at end of file diff --git a/tools/matlab/readrdf.m b/tools/matlab/readrdf.m new file mode 100644 index 0000000000..354ea01f89 --- /dev/null +++ b/tools/matlab/readrdf.m @@ -0,0 +1,142 @@ +function varargout = readrdf(varargin) +% Function to read Radial Distribution Funtion output from LAMMPS +% Input +% 'bin' --> number of bins in rdf histogram +% 'runtime' --> Run length of each of the run commands +% 'step' --> rdf Dump step for each of the run commands +% 'ncol' --> number of columns in the file +% 'final' --> 'YES' indicates that only the average values will be extracted +% If only the averages are needed, you don't need to specify 'runtime', +% 'step' and 'ncol' +% +% Output is in the form of a structure with following variables +% +% if 'final' was given as 'NO' +% .rdf_step_data --> Matrix with RDF for each of the dumped steps +% .rdf_ave_data --> Matrix with average RDF for each RUN +% .rdf_ave_data --> Matrix with average RDF for each RUN +% +% if 'final' was given as 'YES' +% .rdf_ave_data --> Matrix with average RDF for each RUN +% +% Example +% rdf = readrdf('inputfile','bin',100,'runtime',[30000;20000],... +% 'step',[100;100],'ncol',3,'final','no') +% +% 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. + + +rdf_name = varargin{1}; % RDF File name + +% Setting all status checks to zero +bin_status = 0; +runtime_status = 0; +step_status = 0; +ncol_status = 0; +final_status = 'no'; + + + +for id = 2 : 1 : length(varargin) + if strcmpi(varargin{id},'bin') + rdf_bin = varargin{id+1}; % RDF bins + bin_status = 1; + elseif strcmpi(varargin{id},'runtime') + rdf_runtime = varargin{id+1}; % Runtimes + runtime_status = 1; + elseif strcmpi(varargin{id},'step') + rdf_step = varargin{id+1}; % Runtimes + step_status = 1; + elseif strcmpi(varargin{id},'ncol') + ncol = varargin{id+1}; % Runtimes + ncol_status = 1; + elseif strcmpi(varargin{id},'final') + final_status = varargin{id+1}; + end +end + +if ncol_status == 0 + ncol = 3; +end + +% Check for errors in input arguments +if bin_status == 0 + error('No bin specified'); +elseif step_status == 1 && runtime_status == 0 + error('Step size specified without Runtime'); +elseif step_status == 0 && runtime_status == 1 + error('Runtime specified without Step size'); +end +if step_status == 1 && runtime_status == 1 + if length(rdf_runtime) ~= length(rdf_step) + error('Runtime and Step size do not match'); + end +end + +% Preallocating memory if runtime and step size are provided +if step_status == 1 && runtime_status == 1 && strcmpi(final_status,'no') + total_steps = round(sum(rdf_runtime./rdf_step)); + rdf_step_data = zeros(rdf_bin,ncol,total_steps); + rdf_ave_data = zeros(rdf_bin,ncol,length(rdf_runtime)); +elseif strcmpi(final_status,'yes') && step_status == 1 && runtime_status == 1 + rdf_ave_data = zeros(rdf_bin,ncol,length(rdf_runtime)); +end + + +try + rdf_file = fopen(rdf_name,'r'); +catch + error('RDF file not found!'); +end + + + +if strcmpi(final_status,'yes') + run_id = 1; % Run id.. + while feof(rdf_file) ~= 1 + rdf_data = fgetl(rdf_file); + if strcmpi(rdf_data(1:11),'RUN AVERAGE') + fgetl(rdf_file); % to skip the title + for id = 1 : 1 : rdf_bin + rdf_ave_data(id,:,run_id) = str2num(fgetl(rdf_file)); + end + run_id = run_id + 1; + end + + end +else + run_id = 1; + id = 1; + while feof(rdf_file) ~= 1 + rdf_data = fgetl(rdf_file); + if strcmpi(rdf_data(1:8),'TIMESTEP') + timestep(id,1) = str2num(rdf_data(10:length(rdf_data))); + fgetl(rdf_file); % to skip the title + for j = 1 : 1 : rdf_bin + rdf_step_data(j,:,id) = str2num(fgetl(rdf_file)); + end + id = id+1; + elseif strcmpi(rdf_data(1:11),'RUN AVERAGE') + fgetl(rdf_file); % to skip the title + for j = 1 : 1 : rdf_bin + rdf_ave_data(j,:,run_id) = str2num(fgetl(rdf_file)); + end + run_id = run_id + 1; + end + end +end + +fclose(rdf_file); + +if strcmpi(final_status,'no') + out_data.timestep = timestep; + out_data.rdf_step_data = rdf_step_data; + out_data.rdf_ave_data = rdf_ave_data; +else + out_data.rdf_ave_data = rdf_ave_data; +end +varargout{1} = out_data; diff --git a/tools/matlab/scandump.m b/tools/matlab/scandump.m new file mode 100644 index 0000000000..65235561ba --- /dev/null +++ b/tools/matlab/scandump.m @@ -0,0 +1,57 @@ +function [varargout] = scandump(varargin) +% Function to scan LAMMPS dump file +% Input is dumpfile name +% Output is a structure with the following members +% .timestep --> column vector with all the timesteps +% .natoms --> column vector with number of atoms in each timestep +% .position --> column vector with position (for input into readdump) +% .ncol --> column vector with number of columns +% .boxbound --> 3 by 3 by N (number of time steps) 3D array with +% [xlo xhi +% ylo yhi +% zlo zhi]; +% Example +% dump = scandump('dump.LAMMPS'); +% +% See also readdump_one, scandump +% +% 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. + +try + dump = fopen(varargin{1},'r'); +catch + error('Dumpfile not found!'); +end + +i = 1; +out.position(i,1) = 0; + +while ~feof(dump) + id = fgetl(dump); + switch id + case 'ITEM: TIMESTEP' + out.timestep(i,1) = str2num(fgetl(dump)); + case 'ITEM: NUMBER OF ATOMS' + out.Natoms(i,1) = str2num(fgetl(dump)); + case 'ITEM: BOX BOUNDS' + out.boxbound(1,:,i) = str2num(fgetl(dump)); + out.boxbound(2,:,i) = str2num(fgetl(dump)); + out.boxbound(3,:,i) = str2num(fgetl(dump)); + case 'ITEM: ATOMS' + t = str2num(fgetl(dump)); + out.ncol(i,1) = length(t); + for j = 2 : 1 : out.Natoms(i,1) + fgetl(dump); + end + out.position(i+1,1) = ftell(dump); + i = i+1; + end +end + +% Output +varargout{1} = out; +