diff --git a/src/opt_disl.f90 b/src/opt_disl.f90
index 776b3c5..75adcc1 100644
--- a/src/opt_disl.f90
+++ b/src/opt_disl.f90
@@ -530,8 +530,8 @@ module opt_disl
     subroutine vacancy_disloop
         !This subroutine actually creates the dislocation loop.
 
-        real(kind=dp) :: neighbor_dis, temp_box(6), dis, normal_dim
-        integer :: i, j, index, delete_num, delete_index(atom_num)
+        real(kind=dp) :: neighbor_dis, temp_box(6), dis
+        integer :: i, j, index, delete_num, delete_index(atom_num), normal_dim
 
         neighbor_dis = HUGE(1.0_dp)
         index= 0
@@ -560,22 +560,41 @@ module opt_disl
         centroid = r_atom(:,index)
         !Now reset the list for the scanning algorithm
 
+        delete_num = 0
         !Now scan over all atoms again and find the closest vloop_size number of atoms to the initial atom 
-        !that reside on the same plane.
+        !that reside on the same plane. If loop_radius is > 0 then we build a circular vacancy cluster
+        if(loop_radius > 0) then 
+            do i = 1, atom_num
+                !Check to see if it is on the same plane
+                if (in_block_bd(r_atom(:,i), temp_box)) then 
+                    dis = norm2(r_atom(:,i) - centroid)
+                    !Check to see if it is within the loop radius, if so then add it to the delete list
+                    if (dis < loop_radius) then 
+                        delete_num = delete_num + 1
+                        delete_index(delete_num) = i
+                    end if
+                end if
+            end do 
 
-        delete_num = 0
-        do i = 1, atom_num
-            !Check to see if it is on the same plane
-            if (in_block_bd(r_atom(:,i), temp_box)) then 
-                dis = norm2(r_atom(:,i) - centroid)
-                !Check to see if it is within the loop radius, if so then add it to the delete list
-                if (dis < loop_radius) then 
+        !If the loop radius is < 0 we build a square vacancy cluster 
+        else if (loop_radius < 0) then 
+            !Set the new box boundaries
+            do i = 1, 3
+                if (i /= normal_dim) then 
+                    temp_box(2*i) = centroid(i) + abs(loop_radius)
+                    temp_box(2*i-1) = centroid(i) - abs(loop_radius)
+                end if
+            end do 
+
+            !Now delete all atoms with box box boundary which should create a square planar vacancy cluster
+            do i = 1, atom_num
+                !Check to see if the atom is in the bounding cube 
+                if (in_block_bd(r_atom(:,i), temp_box)) then 
                     delete_num = delete_num + 1
                     delete_index(delete_num) = i
                 end if
-            end if
-        end do 
-
+            end do  
+        end if
         !Now delete the atoms
         call delete_atoms(delete_num, delete_index(1:delete_num))