%%1/22/14 - Happy 23rd Birthday me!
%%creats and FCC lattice inside a parallelepiped based on lattice constant, and
%%basis using .xyz notation. 

clc;
clear all;
format long g;
%constants (modify if needed)

%Basis
basis = [0 0 0; 0.5 0.5 0; 0.5 0 0.5; 0 0.5 0.5];

%lattice constant must be defined for temp and pressure, ref papers
latCnst = 1.6255882352;

%Box Dimensions (multiple of lattice constant size, must be integers)

nx = 4;
ny = 4;
nz = 4;

%Atom Types (defaults to 1, change if multiple types and define div
atomType = 1;

% portion of atoms to be devided defaults to 1, for multiple atom types
% define as a fraction, ex. 0.5, 0.25, 0.125 etc. 
div = 1;


%create lattice

%define output name
def = 'in.pos.';
X = 'x';
Y = 'y';
Z = 'z'
x = num2str(nx);
y = num2str(ny);
z = num2str(nz);
 %creat output name in form of box dimensions in.pos.XdimYdimZdim
 
outputname = [def x X y Y z Z];

fid = fopen(outputname,'wt');

%write to file note -1 because matlab is 1:n not 0:n-1
atomNum = 1;
for xdim = 0:nx-1
    for ydim = 0:ny-1
        for zdim = 0:nz-1
            for n = 1:size(basis,1)
                %print atom number
                fprintf(fid, '%d \t', atomNum);
                %print atom type
                
                %print atom positions
                fprintf(fid, '%d \t', atomType);
                fprintf(fid, '%.12f \t', basis(n,1)*latCnst + xdim*latCnst);
                fprintf(fid,'\t');
                fprintf(fid, '%.12f \t',basis(n,2)*latCnst + ydim*latCnst );
                fprintf(fid,'\t');
                fprintf(fid, '%.12f \t',basis(n,3)*latCnst + zdim*latCnst );
                fprintf(fid,'\n');
                atomNum= atomNum+1;
            end
        end
    end
end
fclose(fid);