diff --git a/src/Makefile b/src/Makefile
index 642fce5..87f9306 100644
--- a/src/Makefile
+++ b/src/Makefile
@@ -1,6 +1,6 @@
 FC=ifort
-#FFLAGS=-mcmodel=large -g -O0 -stand f08 -fpe0 -traceback -check bounds,uninit -warn all -implicitnone -no-wrap-margin -heap-arrays
-FFLAGS=-mcmodel=large -Ofast -no-wrap-margin -heap-arrays
+FFLAGS=-mcmodel=large -g -O0 -stand f08 -fpe0 -traceback -check bounds,uninit -warn all -implicitnone -no-wrap-margin -heap-arrays
+#FFLAGS=-mcmodel=large -Ofast -no-wrap-margin -heap-arrays
 MODES=mode_create.o mode_merge.o mode_convert.o
 OPTIONS=opt_disl.o opt_group.o opt_orient.o opt_delete.o 
 OBJECTS=main.o elements.o io.o subroutines.o functions.o atoms.o call_mode.o box.o  $(MODES) $(OPTIONS) call_option.o sorts.o
diff --git a/src/elements.f90 b/src/elements.f90
index 0dffecc..954719f 100644
--- a/src/elements.f90
+++ b/src/elements.f90
@@ -11,7 +11,7 @@ module elements
 
     !Data structures used to represent the CAC elements. Each index represents an element
     character(len=100), allocatable :: type_ele(:) !Element type 
-    integer, allocatable :: size_ele(:), lat_ele(:), sbox_ele(:) !Element size
+    integer, allocatable :: size_ele(:), lat_ele(:), sbox_ele(:), tag_ele(:) !Element size
     real(kind=dp), allocatable :: r_node(:,:,:,:) !Nodal position array
 
     integer, save :: ele_num !Number of elements
@@ -19,7 +19,7 @@ module elements
 
     !Data structure used to represent atoms 
     integer, allocatable ::  type_atom(:)!atom type 
-    integer, allocatable :: sbox_atom(:)
+    integer, allocatable :: sbox_atom(:), tag_atom(:)
     real(kind =dp),allocatable :: r_atom(:,:) !atom position
     integer :: atom_num=0 !Number of atoms
 
@@ -146,7 +146,7 @@ module elements
         
         !Allocate element arrays
         if(n > 0) then 
-            allocate(type_ele(n), size_ele(n), lat_ele(n), sbox_ele(n), r_node(3,max_basisnum, max_ng_node,n), &
+            allocate(type_ele(n), tag_ele(n), size_ele(n), lat_ele(n), sbox_ele(n), r_node(3,max_basisnum, max_ng_node,n), &
                        stat=allostat)
             if(allostat > 0) then 
                 print *, "Error allocating element arrays in elements.f90 because of: ", allostat
@@ -156,7 +156,7 @@ module elements
 
         if(m > 0) then 
             !Allocate atom arrays
-            allocate(type_atom(m), sbox_atom(m), r_atom(3,m), stat=allostat)
+            allocate(type_atom(m), sbox_atom(m), tag_atom(m), r_atom(3,m), stat=allostat)
             if(allostat > 0) then 
                 print *, "Error allocating atom arrays in elements.f90 because of: ", allostat
                 stop
@@ -187,6 +187,11 @@ module elements
             temp_int(ele_size+1:) = 0
             call move_alloc(temp_int, lat_ele)
 
+            allocate(temp_int(n+ele_num+buffer_size))
+            temp_int(1:ele_size) = tag_ele
+            temp_int(ele_size+1:) = 0
+            call move_alloc(temp_int, tag_ele)
+
             allocate(temp_int(n+ele_num+buffer_size))
             temp_int(1:ele_size) = size_ele
             temp_int(ele_size+1:) = 0
@@ -214,6 +219,11 @@ module elements
             temp_int(atom_size+1:) = 0
             call move_alloc(temp_int, type_atom)
 
+            allocate(temp_int(m+atom_num+buffer_size))
+            temp_int(1:atom_size) = tag_atom
+            temp_int(atom_size+1:) = 0
+            call move_alloc(temp_int, tag_atom)
+
             allocate(temp_int(m+atom_num+buffer_size))
             temp_int(1:atom_size) = sbox_atom
             temp_int(atom_size+1:) = 0
@@ -226,15 +236,25 @@ module elements
         end if
     end subroutine
 
-    subroutine add_element(type, size, lat, sbox, r)
+    subroutine add_element(tag, type, size, lat, sbox, r)
         !Subroutine which adds an element to the element arrays
-        integer, intent(in) :: size, lat, sbox
+        integer, intent(in) :: size, lat, sbox, tag
         character(len=100), intent(in) :: type
         real(kind=dp), intent(in) :: r(3, max_basisnum, max_ng_node)
 
+        integer :: newtag
+
         ele_num = ele_num + 1
+
+        if (tag==0) then
+            newtag = ele_num !If we don't assign a tag then pass the tag as the ele_num
+        else
+            newtag = tag
+        end if
+
         !Check to see if we need to grow the arrays
         call grow_ele_arrays(1,0)
+        tag_ele(ele_num) = newtag
         type_ele(ele_num) = type
         size_ele(ele_num) = size
         lat_ele(ele_num) = lat
@@ -245,14 +265,22 @@ module elements
 
     end subroutine add_element
 
-    subroutine add_atom(type, sbox, r)
+    subroutine add_atom(tag, type, sbox, r)
         !Subroutine which adds an atom to the atom arrays
-        integer, intent(in) :: type, sbox
+        integer, intent(in) :: type, sbox, tag
         real(kind=dp), intent(in), dimension(3) :: r
 
+        integer :: newtag
+
         atom_num = atom_num+1
+        if(tag==0) then
+            newtag = atom_num !If we don't assign a tag then pass the tag as the atom_num
+        else
+            newtag = tag
+        end if
         !Check to see if we need to grow the arrays
         call grow_ele_arrays(0,1)
+        tag_atom(atom_num) = tag
         type_atom(atom_num) = type
         r_atom(:,atom_num) = r(:)
         sbox_atom(atom_num) = sbox
diff --git a/src/io.f90 b/src/io.f90
index cd3e3f6..6aafafc 100644
--- a/src/io.f90
+++ b/src/io.f90
@@ -543,13 +543,13 @@ module io
 
         !Write out atoms first
         do i = 1, atom_num
-            write(11,*) i, type_atom(i), sbox_atom(i), r_atom(:,i)
+            write(11,*) tag_atom(i), type_atom(i), sbox_atom(i), r_atom(:,i)
         end do
 
         !Write out the elements, this is written in two stages, one line for the element and then 1 line for 
         !every basis at every node
         do i = 1, ele_num
-            write(11, *) i, lat_ele(i), size_ele(i), sbox_ele(i), type_ele(i)
+            write(11, *) tag_ele(i), lat_ele(i), size_ele(i), sbox_ele(i), type_ele(i)
             do inod = 1, ng_node(lat_ele(i))
                 do ibasis =1, basisnum(lat_ele(i))
                     write(11,*) inod, ibasis, r_node(:, ibasis, inod, i)
@@ -738,19 +738,19 @@ module io
         do i = 1, in_atoms
             read(11,*) j, type, sbox, r(:)
             r = r+newdisplace
-            call add_atom(new_type_to_type(type), sbox+sub_box_num, r)
+            call add_atom(j, new_type_to_type(type), sbox+sub_box_num, r)
         end do
 
         !Read the elements
         do i = 1, in_eles
-            read(11, *) l, type, size, sbox, etype
+            read(11, *) j, type, size, sbox, etype
             do inod = 1, ng_node(type)
                 do ibasis =1, basisnum(type)
-                    read(11,*) j, k, r_innode(:, ibasis, inod)
+                    read(11,*) k, l, r_innode(:, ibasis, inod)
                     r_innode(:,ibasis,inod) = r_innode(:, ibasis, inod) + newdisplace
                 end do
             end do
-            call add_element(etype, size, new_lattice_map(type), sbox+sub_box_num, r_innode)
+            call add_element(j, etype, size, new_lattice_map(type), sbox+sub_box_num, r_innode)
         end do
 
         !Close the file being read
@@ -773,7 +773,7 @@ module io
         real(kind=dp), dimension(3), intent(in) :: displace
         real(kind = dp), dimension(6), intent(out) :: temp_box_bd
 
-        integer :: i, inod, ibasis, j, k, in_eles, in_atoms, ele_types, in_lat_num, in_atom_types, &
+        integer :: i, inod, ibasis, j, k, l, in_eles, in_atoms, ele_types, in_lat_num, in_atom_types, &
                    atom_type_map(100), etype_map(100), etype, lat_type, new_lattice_map(100), &
                    atom_type
         real(kind=dp) :: newdisplace(3), r_in(3,1,8), r_in_atom(3), new_displace(3)
@@ -919,10 +919,10 @@ module io
             do i = 1, in_eles 
                 read(11,*) j, etype, k, lat_type 
                 do inod = 1, 8
-                    read(11, *) j, k, r_in(:,1,inod)
+                    read(11, *) k, l, r_in(:,1,inod)
                     r_in(:,1,inod) = r_in(:,1,inod) + newdisplace
                 end do 
-                call add_element(in_lattype_map(lat_type), etype_map(etype), new_lattice_map(lat_type), sub_box_num + 1, r_in)
+                call add_element(j, in_lattype_map(lat_type), etype_map(etype), new_lattice_map(lat_type), sub_box_num + 1, r_in)
             end do
         end if
 
@@ -937,7 +937,7 @@ module io
             do i = 1, in_atoms
                 read(11,*) j, k, atom_type, r_in_atom(:)
                 r_in_atom = r_in_atom + newdisplace
-                call add_atom(atom_type_map(atom_type), sub_box_num + 1, r_in_atom)
+                call add_atom(j,atom_type_map(atom_type), sub_box_num + 1, r_in_atom)
             end do 
             !Close file
             close(11)
diff --git a/src/mode_create.f90 b/src/mode_create.f90
index 5c347ab..454fcb7 100644
--- a/src/mode_create.f90
+++ b/src/mode_create.f90
@@ -109,7 +109,7 @@ module mode_create
                 box_bd(2*i) = maxval(r_node_temp(i,:,:))+10.0_dp**-6.0_dp
                 box_bd(2*i-1) = minval(r_node_temp(i,:,:)) - 10.0_dp**-6.0_dp
             end do
-            call add_element(element_type, esize, 1, 1, r_node_temp)
+            call add_element(0,element_type, esize, 1, 1, r_node_temp)
         end if
         
         !If we passed the dup_flag or dim_flag then we have to convert the lattice points and add them to the atom/element arrays
@@ -136,7 +136,7 @@ module mode_create
             if(lat_atom_num > 0) then 
                 do i = 1, lat_atom_num
                     do ibasis = 1, basisnum(1)
-                        call add_atom(basis_type(ibasis, 1), 1, (r_atom_lat(:,i)*lattice_parameter)+basis_pos(:,ibasis))
+                        call add_atom(0,basis_type(ibasis, 1), 1, (r_atom_lat(:,i)*lattice_parameter)+basis_pos(:,ibasis))
                     end do
                 end do
                 deallocate(r_atom_lat)
@@ -150,7 +150,7 @@ module mode_create
                         end do
                     end do 
 
-                    call add_element(element_type, elat(i), 1, 1, r_node_temp)
+                    call add_element(0,element_type, elat(i), 1, 1, r_node_temp)
                 end do 
             end if
         end if
diff --git a/src/opt_group.f90 b/src/opt_group.f90
index 281ab28..c24b857 100644
--- a/src/opt_group.f90
+++ b/src/opt_group.f90
@@ -564,7 +564,7 @@ module opt_group
                 !here as well to make sure they are in the box
                 do j = 1, basisnum(lat_ele(ie))*size_ele(ie)**3
                     call apply_periodic(r_interp(:,j))
-                    call add_atom(type_interp(j), sbox_ele(ie), r_interp(:,j))
+                    call add_atom(0,type_interp(j), sbox_ele(ie), r_interp(:,j))
                 end do 
             end do
             !Once all atoms are added we delete all of the elements
@@ -800,7 +800,7 @@ module opt_group
                                         !Add the element, for the sbox we just set it to the same sbox that we get the orientation 
                                         !from. In this case it is from the sbox of the first atom in the group.
                                         new_ele = new_ele+1
-                                        call add_element(remesh_ele_type, working_esize, ilat, sbox_atom(atom_index(1)),r_new_node)
+                                        call add_element(0,remesh_ele_type, working_esize, ilat, sbox_atom(atom_index(1)),r_new_node)
 
                                     end if
                                 end if
@@ -821,7 +821,7 @@ module opt_group
                                 lat_points(ix,iy,iz) = .false.
                                r = matmul(orient, matmul(fcc_mat, vlat))*lapa(ilat)
                                new_atom=new_atom+1
-                               call add_atom(basis_type(1,ilat), is, r)
+                               call add_atom(0,basis_type(1,ilat), is, r)
                             end if
                         end do
                     end do