% This program is designed to give atomic postions of an FCC material.

% The intent is to model a bicrystal by simulating a twist grain boundary

% along the 110 plane



% Experimentally Determined Atomic Radii for common FCC materials

% Al=1.431, Copper=1.278, Gold=1.442, Lead=1.750, Nickel=1.246

% Platinum=1.387, Silver=1.445;  All units are in Angstroms = 10^-10 m

% Theoretical calculation of the lattice spacing a from the atomic radii:

% a=r*2^1.5;

%

clear

% Dialog box asks user to enter parameters

prompts =  {'Enter number of repeating stacking planes for length of the PBC:',

    'Enter angle of twist in degrees (90, 43.6, 36.9, 31.9, 28.0, 22.7, 18.9, 16.3, 12.7, 10.4, 8.8, 0):',

    'Enter lattice spacing (in Angstroms) for material of interest:',

    'Enter file name for output of the XYZ coordinate position file:'};

defaults = {'5','38.9','3.524','FCC_389_5.txt'};

answer = inputdlg(prompts,'Enter inputs for Common Site Lattice Algorithm',1,defaults);

L              = str2num(answer{1});

phi            = str2num(answer{2});

a              = str2num(answer{3}); 

filename       = answer{4};

%

% Outputs the various X-Y box size according to the twist angle 

% as governed by the Pathagorean triplets

if phi==90

    N=5;

end

if phi==43.6

    N=29;

end

if phi==36.9

    N=5;

end

if phi==31.9

    N=53;

end

if phi==28.0

    N=17;

end

if phi==22.7

    N=13;

end

if phi==18.9

    N=37;

end

if phi==16.3

    N=25;

end

if phi==12.7

    N=41;

end

if phi==10.4

    N=61;

end

if phi==8.8

    N=85;

end

if phi==0

    N=5;

end

if phi==38.9

    N=5;

end

if phi==54.7

    N=5;

end


if phi==70.5

    N=20;

end

% Error associated with removing atoms given their atomic radii

err=0.1;

%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Create the atomic positions for the untwisted lattice

% This lattice is much bigger than needed

% And the extra atoms will be removed after superimposing

% the twisted lattice on top

C1=zeros((L)*(2*N+1)*(4*N+2)+2*N*L*(4*N+1),3);

% Populate the first grain or lattice structure

m=1;

f=sqrt(2)/2;

for k=1:L

    for j=1:2*N+1

        for i=1:4*N+2

            C1(m,1)=(i-1)*a*f-(2*a*f*(N-1.5));

            C1(m,2)=(j-1)*a;

            C1(m,3)=(k-1)*a*f;

            m=m+1;

        end

    end

end

i=1;

j=1;

k=1;

for k=1:L

    for j=1:2*N

        for i=1:4*N+1

            C1(m,1)=(i-0.5)*a*f-(2*a*f*(N-1.5));

            C1(m,2)=(j-0.5)*a;

            C1(m,3)=(k-0.5)*a*f;

            m=m+1;

        end

    end

end

%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Create the atomic positions for the twisted lattice

C2=zeros((L)*(N+1)*(2*N+1)+L*N*2*N,3);

% Populate the first grain or lattice structure

i=1;

j=1;

k=1;

m=1;

for k=1:L

    for j=1:N+1

        for i=1:2*N+1

            C2(m,1)=(i-1)*a*f;

            C2(m,2)=(j-1)*a;

            C2(m,3)=(k-1)*a*f;

            m=m+1;

        end

    end

end

i=1;

j=1;

k=1;

for k=1:L

    for j=1:N

        for i=1:2*N

            C2(m,1)=(i-0.5)*a*f;

            C2(m,2)=(j-0.5)*a;

            C2(m,3)=(k-0.5)*a*f;

            m=m+1;

        end

    end

end

C2max=max(size(C2));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Twist the second lattice structure

% Define the Rotation Matrix and convert phi to radians

phi_r=pi*phi/180;

R=[cos(phi_r) sin(phi_r) 0; -sin(phi_r) cos(phi_r) 0; 0 0 1];

m=1;

% Rotates or twists the structure and shifts it's Z component downward

for m=1:C2max

    C2(m,:)=C2(m,:)*R;

    C2(m,3)=C2(m,3)-(L)*a*f;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Create the atomic positions for the twisted lattice

C4=zeros((L+1)*(N+1)*(2*N+1)+L*N*2*N,3);

% Populate the first grain or lattice structure

i=1;

j=1;

k=1;

m=1;

for k=1:L+1

    for j=1:N+1

        for i=1:2*N+1

            C4(m,1)=(i-1)*a*f;

            C4(m,2)=(j-1)*a;

            C4(m,3)=(k-1)*a*f;

            m=m+1;

        end

    end

end

i=1;

j=1;

k=1;

for k=1:L

    for j=1:N

        for i=1:2*N

            C4(m,1)=(i-0.5)*a*f;

            C4(m,2)=(j-0.5)*a;

            C4(m,3)=(k-0.5)*a*f;

            m=m+1;

        end

    end

end

C4max=max(size(C4));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Twist the third lattice structure


m=1;

% Rotates or twists the structure and shifts it's Z component downward

for m=1:C4max

    C4(m,:)=C4(m,:)*R;

    C4(m,3)=C4(m,3)+(L)*a*f;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Removes all atoms from the first lattice that are not within

% the periodic boundary box size corresponding to the second twisted lattice

j=1;

i=1;

Num_atoms=max(size(C1));

C3=zeros(1,3);

%

% Loop over all the atoms and if they are 'in bounce', 

% move them to array C3 

err2=err;

for j=1:Num_atoms

    if j>Num_atoms

        break

    end

    % Account for two special cases to avoid division by 0

    if phi==0

        SINE_phi_r=1;

        COSINE_phi_r=cos(phi_r);

        X1=0;

        X2=C1(j,1);

        Y1=C1(j,2);

        Y2=0;

    elseif phi==90

        SINE_phi_r=sin(phi_r);

        COSINE_phi_r=1;

        X1=-C1(j,1);

        X2=0;

        Y1=0;

        Y2=C1(j,2);

    else

        SINE_phi_r=sin(phi_r);

        COSINE_phi_r=cos(phi_r);

        X1=-tan(phi_r)*C1(j,1);

        X2=((tan(phi_r))^-1)*C1(j,1);

        Y1=C1(j,2);

        Y2=C1(j,2);

    end

    Cond1=Y1+X1;

    Cond2=Y2+X2;

    Line1=0;

    Line3=((2*N)*a*f)/SINE_phi_r;

    Line4=((N)*a)/COSINE_phi_r;

    if Cond2>=Line1

        if Cond1>=Line1

            if Cond2<=Line3

                if Cond1<=Line4

                    C3(i,:)=C1(j,:);

                    i=i+1;

                end

            end

        end

    end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Concatinates the first and second lattice into one array of XYZ coords

C=vertcat(C3,C2,C4);

Cmax=max(size(C));

i=1;

% Rotates the entire array, simulating the grain boundary, so that

% it is aligned with the original coord axis and in the first quadrant

for i=1:Cmax

    C(i,:)=C(i,:)*R^-1;

    C(i,3)=C(i,3)+L*a*f;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Determines the 8 corners of the PBC box in order to apply boundary

% conditions in the MD simulation

Xmin=min(C(:,1));

Xmax=max(C(:,1));

Ymin=min(C(:,2));

Ymax=max(C(:,2));

Zmin=min(C(:,3));

Zmax=max(C(:,3));

%err1=err+a*f;

err1=err;

i=1;

for i=1:Cmax

    if C(i,:)>[Xmin-err1 Ymin-err1 Zmax-err] & C(i,:)<[Xmin+err1 Ymin+err1 Zmax+err]

        Atom1=i

    end

     if C(i,:)>[Xmax-2*err1 Ymin-err1 Zmax-err] & C(i,:)<[Xmax+2*err1 Ymin+err1 Zmax+err]

        Atom2=i

    end   

    if C(i,:)>[Xmax-2*err1 Ymax-2*err1 Zmax-err] & C(i,:)<[Xmax+2*err1 Ymax+2*err1 Zmax+err]

        Atom3=i

    end

        if C(i,:)>[Xmin-err1 Ymax-2*err1 Zmax-err] & C(i,:)<[Xmin+err1 Ymax+2*err1 Zmax+err]

        Atom4=i

    end

        if C(i,:)>([Xmin Ymin Zmin]-err) & C(i,:)<([Xmin Ymin Zmin]+err)

        Atom5=i

    end

        if C(i,:)>[Xmax-2*err1 Ymin-err1 Zmin-err] & C(i,:)<[Xmax+2*err1 Ymin+err1 Zmin+err]

        Atom6=i

    end

        if C(i,:)>[Xmax-2*err1 Ymax-2*err1 Zmin-err] & C(i,:)<[Xmax+2*err1 Ymax+2*err1 Zmin+err]

        Atom7=i

    end

        if C(i,:)>[Xmin-err1 Ymax-2*err1 Zmin-err] & C(i,:)<[Xmin+err1 Ymax+2*err1 Zmin+err]

        Atom8=i

    end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Plots the response: a 3D graph and 3-2D graphs 

plot3(C(:,1),C(:,2),C(:,3),'.')

figure

plot(C(:,1),C(:,2),'.')

figure

plot(C(:,1),C(:,3),'.')

figure

plot(C(:,2),C(:,3),'.')

i=1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Saves the atoms' XYZ coord pairs to a file of the user's choice

fid=fopen(filename,'w');

for i=1:Cmax

    fprintf(fid, '%6.4f  %6.4f  %6.4f \n',C(i,1),C(i,2),C(i,3));

end

fclose(fid);

Zmax

Number_of_Atoms=Cmax