diff --git a/AddRectAtomicArray.m b/AddRectAtomicArray.m
new file mode 100644
index 0000000..93de8d1
--- /dev/null
+++ b/AddRectAtomicArray.m
@@ -0,0 +1,55 @@
+function AddRectAtomicArray(LAtoms, WAtoms, X0, Y0, VX0, VY0, InitDist, Temp, Type)
+global C
+global x y AtomSpacing
+global nAtoms
+global AtomType Vx Vy Mass0 Mass1
+
+if Type == 0            %sets mas based on the type of atom set in MD.m (global variable)
+    Mass = Mass0;
+else
+    Mass = Mass1;
+end
+
+%sets the lenght and with based on the number of atoms provided when
+%calling this function
+L = (LAtoms - 1) * AtomSpacing;
+W = (WAtoms - 1) * AtomSpacing;
+
+numAtoms = LAtoms * WAtoms;
+
+xp(1, :) = linspace(0, L, LAtoms);
+yp(1, :) = linspace(0, W, WAtoms);
+
+x(nAtoms + 1:nAtoms+LAtoms) = xp-L/2;
+y(nAtoms + 1:nAtoms+LAtoms) = yp(1)-W/2;
+
+for i = 1:WAtoms-1
+    x(nAtoms + i * LAtoms + 1:nAtoms + (i + 1) * LAtoms) = xp - L / 2;
+    y(nAtoms + i * LAtoms + 1:nAtoms + (i + 1) * LAtoms) = yp(i + 1) - W / 2;
+end
+
+x(nAtoms + 1:nAtoms + numAtoms) = x(nAtoms + 1:nAtoms + numAtoms) + ...
+    (rand(1, numAtoms) - 0.5) * AtomSpacing * InitDist + X0;
+y(nAtoms + 1:nAtoms + numAtoms) = y(nAtoms + 1:nAtoms + numAtoms) + ...
+    (rand(1, numAtoms) - 0.5) * AtomSpacing * InitDist + Y0;
+
+AtomType(nAtoms + 1:nAtoms + numAtoms) = Type;
+
+if Temp == 0
+    Vx(nAtoms + 1:nAtoms + numAtoms) = 0;
+    Vy(nAtoms + 1:nAtoms + numAtoms) = 0;
+else
+    std0 = sqrt(C.kb * Temp / Mass);
+
+    Vx(nAtoms + 1:nAtoms + numAtoms) = std0 * randn(1, numAtoms);
+    Vy(nAtoms + 1:nAtoms + numAtoms) = std0 * randn(1, numAtoms);
+end
+
+Vx(nAtoms + 1:nAtoms + numAtoms) = Vx(nAtoms + 1:nAtoms + numAtoms) - ...
+    mean(Vx(nAtoms + 1:nAtoms + numAtoms)) + VX0;
+Vy(nAtoms + 1:nAtoms + numAtoms) = Vy(nAtoms + 1:nAtoms + numAtoms) - ...
+    mean(Vy(nAtoms + 1:nAtoms + numAtoms)) + VY0;
+
+nAtoms = nAtoms + numAtoms;
+
+end
diff --git a/GetForces.m b/GetForces.m
new file mode 100644
index 0000000..edc7916
--- /dev/null
+++ b/GetForces.m
@@ -0,0 +1,60 @@
+function GetForces(PhiCutoff,Epsilon,sigma)
+global x y nAtoms Fx Fy Phi
+global step applyForce LAtoms WAtoms
+
+
+squish = false;             %%initially squish is set to false
+if applyForce               %%if apply force is true, force will initally be applied - tried to see how interactions are after large force is applied inwards
+    if step == 0            %%only applies the force initially (time 0)
+        squish = true;             %%will apply force fo 1e-9 to top and bottom electrons (facing inwards) -  atomic array will fall apart
+    end
+end
+
+% n = 1;
+for i = 1:nAtoms            %%cycles through for each atoms
+    Fx(i) = 0;              %%initially set everything for that atom to zero
+    Fy(i) = 0;
+    Phi(i) = 0;
+%     if i == 10
+%         plot(x(i),y(i),'o','markers',48);
+%         hold on
+%     end
+
+    for j = 1:nAtoms            %cycles through atoms to calcuate force on the number i atoms
+        if i == j, continue; end        %%if i = j it will pass to the next iteration of the for loop - can't ahve a force from itself
+        dx = x(i) - x(j);
+        dy = y(i) - y(j);
+        r = sqrt(dx^2 + dy^2);
+
+        if r > PhiCutoff, continue, end         %if atom is too far away, it does not have an effect on the $i atom
+
+        [aPhi dPhidr] = LJPot(r, Epsilon, sigma);
+
+        ang = atan2(dy, dx);
+%         dx =  r*cos(ang)
+%         dy =  r*sin(ang)
+
+        dFx = - dPhidr * cos(ang);
+        dFy = - dPhidr * sin(ang);
+
+%         if i == 10
+%             plot(x(j),y(j),'ro','markers',48);
+%         end
+
+        Phi(i) = Phi(i) + aPhi;
+        Fx(i) = Fx(i) + dFx;
+        Fy(i) = Fy(i) + dFy;   
+       
+        if squish                   %%apply force to top and bottom layers (force inwards) - ignores the calculated forces
+            if i > (LAtoms*WAtoms - LAtoms)
+                Fy(i) = -1e-9; 
+            elseif i < (LAtoms+1)
+                Fy(i) = 1e-9;
+            end
+    end
+end
+
+hold off
+
+end
+
diff --git a/InitBlock.m b/InitBlock.m
new file mode 100644
index 0000000..4dcfdc5
--- /dev/null
+++ b/InitBlock.m
@@ -0,0 +1,36 @@
+global LAtoms WAtoms
+
+%define globals to set number of atoms in rectangle
+LAtoms = 10;        
+WAtoms = 10;
+
+doPlot = 1;
+dt = 15e-15;
+TStop = 1000 * dt;
+InitDist = 0.0;
+Method = 'VE'; % VE -- verlot; FD -- Forward Difference
+
+Mass0 = 14 * C.am; % Silicon
+Mass1 = 5 * C.am; % Argon
+
+AtomSpacing = 0.5430710e-9;
+LJSigma = AtomSpacing / 2^(1/6);
+LJEpsilon = 1e-21;
+
+PhiCutoff = 3 * AtomSpacing * 1.1;      %cut-off set to more than three times the atomic spacing - these atoms will not interact with one another
+
+T = 30;                                 %temperature set to 30 degrees
+
+AddRectAtomicArray(LAtoms, WAtoms, 0, 0, 0, 0, 0, T, 0);
+
+%sets up the plot sizing
+Size = 8 * AtomSpacing;
+Limits = [-Size +Size -Size +Size]; % square is good
+PlDelt = 5*dt;
+MarkerSize = 10;
+PlotFile = 'Block.gif';
+doPlotImage = 1;
+PlotSize = [100, 100, 1049, 1049];
+
+ScaleV = .2e-11;
+ScaleF = 20;
diff --git a/MD.m b/MD.m
new file mode 100644
index 0000000..975632e
--- /dev/null
+++ b/MD.m
@@ -0,0 +1,201 @@
+% clear all
+clearvars
+clearvars -GLOBAL
+close all
+format shorte
+
+%define global variables used
+set(0, 'DefaultFigureWindowStyle', 'docked')
+global C
+global Vx Vy x y Fx Fy AtomSpacing
+global Phi nAtoms time Mass0 Mass1 Pty0in Pty1in
+global LJEpsilon LJSigma Phi0 AtomType
+global MinX MaxX MinY MaxY PhiTot KETot
+global nAtoms0 nAtoms1 T T0 T1 MarkerSize
+global doPlotImage PlotCount map im PlotSize ScaleV ScaleF
+global step applyForce LAtoms WAtoms
+
+%set globals
+C.q_0 = 1.60217653e-19;             % electron charge
+C.hb = 1.054571596e-34;             % Dirac constant
+C.h = C.hb * 2 * pi;                    % Planck constant
+C.m_0 = 9.10938215e-31;             % electron mass
+C.kb = 1.3806504e-23;               % Boltzmann constant
+C.eps_0 = 8.854187817e-12;          % vacuum permittivity
+C.mu_0 = 1.2566370614e-6;           % vacuum permeability
+C.c = 299792458;                    % speed of light
+C.g = 9.80665; %metres (32.1740 ft) per sē
+C.am = 1.66053892e-27;
+
+step = 0;                   %used to count the step number (how many times getForce is calculated)
+%applyForce = true;          %%true - apply force to top and bottom layers at step 0 in the getForce function -works with InitBlock only (any size rectangle)
+applyForce = false;             %%false - no force is applied
+
+MaxX = 0;
+MinX = 0;
+MaxY = 0;
+MinY = 0;
+nAtoms = 0;
+MarkerSize = 12;
+Limits = [];
+doPlot = 1;
+doPlotImage = 0;
+PlotCount  = 0;
+PlotFile = 'image.gif';
+PlotSize = [100, 100, 1049, 895];
+ScaleV = 0;
+ScaleF = 0;
+
+% Simulation initiallization - varies for different geometries
+ InitBlock
+% InitBlock0
+% InitBlock0FD
+% InitVStream
+% InitHCP
+% InitHCPBlob
+% InitVStreamHCP
+
+MaxX = max(x) * 1.5;
+MinX = min(x) * 1.5;
+
+MaxY = max(y) * 1.5;
+MinY = min(y) * 1.5;
+%define arrays/vectors used - each 'point' in vector represents an atom
+Fx = zeros(1, nAtoms);
+Fy = zeros(1, nAtoms);
+Phi = zeros(1, nAtoms);
+dx = zeros(1, nAtoms);
+dy = zeros(1, nAtoms);
+dvx = zeros(1, nAtoms);
+dvy = zeros(1, nAtoms);
+
+Pty0in = AtomType == 0;
+Pty1in = AtomType == 1;
+
+nAtoms0 = sum(Pty0in);
+nAtoms1 = sum(Pty1in);
+
+%find inital forces on atoms
+GetForces(PhiCutoff,LJEpsilon,LJSigma);
+Phi0 = sum(Phi) / nAtoms / 2;
+
+V2 = Vx.*Vx + Vy.*Vy;
+% KEc = 1/2*Mass*mean(V2);
+% Tc = KEc/C.kb;
+
+t = 0;
+c = 1;
+time(c) = 0;
+
+PhiTot(c) = sum(Phi) / 2;
+
+V2_0 = (Vx(Pty0in).*Vx(Pty0in) + Vy(Pty0in).*Vy(Pty0in));
+if nAtoms1
+    V2_1 = (Vx(Pty1in).*Vx(Pty1in) + Vy(Pty1in).*Vy(Pty1in));
+else
+    V2_1 = 0;
+end
+
+KE0 = mean(V2_0) * Mass0 * 0.5;
+KE1 = mean(V2_1) * Mass1 * 0.5;
+
+KETot(c) = (KE0 * nAtoms0 + KE1 * nAtoms1);
+T(c) = KETot(c) / nAtoms / C.kb;
+T0(c) = KE0 / C.kb;
+T1(c) = KE1 / C.kb;
+
+PlotVars(c, Limits);
+
+xp = x - dt * Vx;
+xpp = x - 2 * dt * Vx;
+yp = y - dt * Vy;
+ypp = y - 2 * dt * Vy;
+
+Plt0 = PlDelt;
+
+while t < TStop         %cycles through a set number of steps
+
+    %     F = ma
+    %     F = m dv/dt
+
+    GetForces(PhiCutoff,LJEpsilon,LJSigma);
+    step = step + 1;            %%keeps track of the number of steps
+    % Forward difference
+    if Method == 'FD'
+        %     dv = F/m dt
+        %     x = Vx * dt + F/m (dt)^2 / 2
+
+        dvx(Pty0in) = Fx(Pty0in) * dt / Mass0;
+        dvx(Pty1in) = Fx(Pty1in) * dt / Mass1;
+
+        Vx = Vx + dvx;
+        dx(Pty0in) = Vx(Pty0in) * dt + Fx(Pty0in) * dt^2 / 2 / Mass0;
+        dx(Pty1in) = Vx(Pty1in) * dt + Fx(Pty1in) * dt^2 / 2 / Mass1;
+
+        dvy(Pty0in) = Fy(Pty0in) * dt/Mass0;
+        dvy(Pty1in) = Fy(Pty1in) * dt/Mass1;
+
+        Vy = Vy + dvy;
+
+        dy(Pty0in) = Vy(Pty0in) * dt + Fy(Pty0in) * dt^2 / 2 / Mass0;
+        dy(Pty1in) = Vy(Pty1in) * dt + Fy(Pty1in) * dt^2 / 2 / Mass1;
+
+        x = xp + dx;
+        y = yp + dy;
+
+    elseif Method == 'VE'
+
+        x(Pty0in) = -xpp(Pty0in) + 2 * xp(Pty0in) + dt^2 / Mass0 * Fx(Pty0in);
+        x(Pty1in) = -xpp(Pty1in) + 2 * xp(Pty1in) + dt^2 / Mass1 * Fx(Pty1in);
+
+        y(Pty0in) = -ypp(Pty0in) + 2 * yp(Pty0in) + dt^2 / Mass0 * Fy(Pty0in);
+        y(Pty1in) = -ypp(Pty1in) + 2 * yp(Pty1in) + dt^2 / Mass1 * Fy(Pty1in);
+
+        Vx = (x - xpp) / (2 * dt);
+        Vy = (y - ypp) / (2 * dt);
+    end
+
+    xpp = xp;
+    ypp = yp;
+
+    xp = x;
+    yp = y;
+
+    c = c + 1;
+    t  = t + dt;
+    time(c) = t;
+
+
+    PhiTot(c) = sum(Phi)/2;
+    V2_0 = (Vx(Pty0in).*Vx(Pty0in)+Vy(Pty0in).*Vy(Pty0in));
+    if nAtoms1
+        V2_1 = (Vx(Pty1in).*Vx(Pty1in)+Vy(Pty1in).*Vy(Pty1in));
+    else
+        V2_1 = 0;
+    end
+
+    KE0 = mean(V2_0) * Mass0 * 0.5;
+    KE1 = mean(V2_1) * Mass1 * 0.5;
+
+    KETot(c) = (KE0 * nAtoms0 + KE1 * nAtoms1);
+    T(c) = KETot(c) / nAtoms / C.kb;
+    T0(c) = KE0 / C.kb;
+    T1(c) = KE1 / C.kb;
+
+    if t > Plt0
+        fprintf('time: %g (%5.2g %%)\n', t, t / TStop * 100);
+
+        PlotVars(c, Limits);
+
+        Plt0 = Plt0 + PlDelt;
+        pause(0.000001)
+    end
+
+end
+
+
+if doPlotImage
+    imwrite(im, map, PlotFile, 'DelayTime', 0.05, 'LoopCount', inf);
+end
+
+