diff --git a/src/io.f90 b/src/io.f90
index f929e29..1e1d956 100644
--- a/src/io.f90
+++ b/src/io.f90
@@ -161,7 +161,7 @@ module io
         !Calculate total atom number
         write_num = atom_num
         do i = 1,ele_num
-            if(type_ele(i) == 'fcc') write_num = write_num + size_ele(i)**3
+        if((type_ele(i) == 'fcc').or.(type_ele(i) == 'bcc')) write_num = write_num + size_ele(i)**3
         end do
         !Write total number of atoms + elements
         write(11, '(i16, a)') write_num, ' atoms'
@@ -196,7 +196,7 @@ module io
         do i = 1, ele_num
             call interpolate_atoms(type_ele(i), size_ele(i), lat_ele(i), r_node(:,:,:,i), type_interp, r_interp)
             select case(trim(adjustl(type_ele(i))))
-            case('fcc')
+            case('fcc','bcc')
                 do iatom = 1, basisnum(lat_ele(i))*size_ele(i)**3
                     interp_num = interp_num+1
                     call apply_periodic(r_interp(:,iatom))
@@ -949,4 +949,66 @@ module io
             call set_max_esize
         end if
     end subroutine read_pycac
+
+    subroutine read_cac(file, displace, temp_box_bd)
+        !This subroutine is used to read .cac files which are used with the lammpsCAC format
+        character(len=100), intent(in) :: file
+        real(kind=dp), dimension(3), intent(in) :: displace
+        real(kind = dp), dimension(6), intent(out) :: temp_box_bd
+        
+        character(len=100) :: textholder, element_type, esize
+        character(len=2) :: atom_species
+       integer :: i, j, ele_in, type_in, type_map(10), in_basis, node_types(10,8), inod, ibasis, in_basis_types(10)
+                
+        real(kind=dp) :: mass, r_in(3,10,8)
+        
+
+        !First open the file
+        open(unit=11, file=trim(adjustl(file)), action='read',position='rewind')
+
+        !Read header information
+        read(11, *) textholder
+        read(11, *) textholder
+
+        !Read number of elements
+        read(11, *) ele_in, textholder 
+        read(11, *) type_in, textholder
+
+        !Read box_boundaries
+        read(11,*) textholder
+        read(11,*) box_bd(1:2), texholder
+        read(11,*) box_bd(3:4), texholder
+        read(11,*) box_bd(5:6), texholder
+
+        !Read useless information
+        read(11,*) textholder 
+        read(11,*) textholder 
+
+        !Read atomic masses
+        do i = 1, type_in
+            read(11,*) j, mass, textholder
+            call ATOMMASSSPECIES(mass, atom_species) 
+            call add_atom_type(atom_species, type_map(i))
+        end do
+
+        !Read useless info
+        do i = 1, 3
+            read(11,*) textholder
+        end do
+
+        !Start the reading loop
+        do i = 1, ele_in
+            read(11,*) j, ele, element_type, in_basis, esize
+            select case(trim(adjustl(element_type)))
+            case('Eight_Node')
+                    !Read in all the data
+                    do j = 1, 8*in_basis 
+                        read(11, *) inod, ibasis, in_basis_types(ibasis), r_in(:,ibasis,inod)
+                    end do
+
+                    !
+            end select
+        end do
+
+    end subroutine read_cac
 end module io