From 54aa50b605f5e8735f944cbf71e9108505417bee Mon Sep 17 00:00:00 2001
From: Alex Selimov <aselimov3@gatech.edu>
Date: Fri, 29 Nov 2019 13:36:19 -0500
Subject: [PATCH 01/17] Working writing to lammpsCAC format

---
 src/io.f90          | 161 +++++++++++++++++++++++++++++++++++++++-----
 src/mode_create.f90 |   2 -
 2 files changed, 146 insertions(+), 17 deletions(-)

diff --git a/src/io.f90 b/src/io.f90
index 7850ee2..948fd52 100644
--- a/src/io.f90
+++ b/src/io.f90
@@ -6,8 +6,8 @@ module io
 
     implicit none
 
-    integer :: outfilenum = 0
-    character(len=100) :: outfiles(10)
+    integer :: outfilenum = 0, infilenum = 0
+    character(len=100) :: outfiles(10), infiles(10)
 
    public 
    contains
@@ -41,6 +41,7 @@ module io
                 if((scan(overwrite, "n") > 0).or.(scan(overwrite, "N") > 0)) then 
                     print *, "Please specify a new filename with extension:"
                     read(*,*) temp_outfile
+                    cycle
                 else if((scan(overwrite, "y") > 0).or.(scan(overwrite, "Y") > 0)) then 
                     continue
                 else
@@ -56,21 +57,13 @@ module io
                 cycle
             end if
             select case(temp_outfile(scan(temp_outfile,'.',.true.)+1:))
-            case('xyz')
+            case('xyz','lmp','vtk','cac')
                 outfilenum=outfilenum+1
                 outfiles(outfilenum) =  temp_outfile
                 exit
-            case('lmp')
-                outfilenum=outfilenum+1
-                outfiles(outfilenum) = temp_outfile
-                exit
-            case('vtk')
-                outfilenum=outfilenum+1
-                outfiles(outfilenum)=temp_outfile
-                exit
             case default    
-                print *, "File type: ", trim(temp_outfile(scan(temp_outfile,'.',.true.):)), "not currently accepted. ", &
-                         "please input a filename with extension from following list: xyz, lmp, vtk."
+                print *, "File type: ", trim(temp_outfile(scan(temp_outfile,'.',.true.):)), " not currently accepted. ", &
+                         "please input a filename with extension from following list: xyz, lmp, vtk, cac."
                 read(*,*) temp_outfile
 
             end select
@@ -78,6 +71,52 @@ module io
 
     end subroutine get_out_file 
 
+    subroutine get_in_file(filename)
+
+        implicit none
+
+        character(len=100), intent(in) :: filename
+        character(len=100) :: temp_infile
+        logical :: file_exists
+
+        !If no filename is provided then this function is called with none and prompts user input
+        if (filename=='none') then 
+            print *, "Please specify a filename with extension to read in:"
+            read(*,*) temp_infile
+        else
+            temp_infile = filename
+        end if
+
+        !Infinite loop which only exists if user provides valid filetype
+        do while(.true.)
+
+            !Check to see if file exists, if it doesn't then ask the user for another input
+            inquire(file=trim(temp_infile), exist=file_exists)
+            if (.not.file_exists) then 
+                print *, "The file ", temp_infile, "does not exist, please input an existing file to read in."
+                read(*,*) temp_infile
+                cycle
+            end if
+
+            if (scan(temp_outfile,'.',.true.) == 0) then 
+                print *, "No extension included on filename, please type a full filename that includes an extension."
+                read(*,*) temp_infile
+                cycle
+            end if
+            select case(temp_outfile(scan(temp_outfile,'.',.true.)+1:))
+            case('cac')
+                infilenum=infilenum+1
+                infiles(infilenum) =  temp_infile
+                exit
+            case default    
+                print *, "File type: ", trim(temp_infile(scan(temp_outfile,'.',.true.):)), " not currently accepted. ", &
+                         "please input a filename with extension from following list: cac."
+                read(*,*) temp_infile
+
+            end select
+        end do
+
+    end subroutine get_in_file 
 
     subroutine write_out
         !This subroutine loops over alll of the outfile types defined and calls the correct writing subroutine
@@ -96,6 +135,8 @@ module io
                 call write_lmp(outfiles(i))
             case('vtk')
                 call write_vtk(outfiles(i))
+            case('cac')
+                call write_lmpcac(outfiles(i))
             case default
                 print *, "The extension ", trim(adjustl(outfiles(i)(scan(outfiles(i),'.',.true.)+1:))), &
                          " is not accepted for writing. Please select from: xyz and try again"
@@ -201,6 +242,71 @@ module io
         end do
     end subroutine write_lmp
 
+    subroutine write_lmpcac(file)
+        !This subroutine writes out a .lmp style dump file
+        character(len=100), intent(in) :: file
+        integer :: write_num, i, inod, ibasis
+        real(kind=dp) :: mass
+
+1       format(i16, '   Eight_Node', 4i16)
+2       format(i16, '   Atom', 4i16)
+3       format(3i16,3f23.15)
+
+        open(unit=11, file=trim(adjustl(file)), action='write', status='replace',position='rewind')
+
+        !Comment line
+        write(11, '(a)') '# CAC input file made with cacmb'
+        write(11, '(a)')
+        !Calculate total atom number
+        write_num = atom_num + ele_num
+
+        !Write total number of atoms + elements
+        write(11, '(i16, a)') write_num, ' cac elements'
+        !Write number of atom types 
+        write(11, '(i16, a)') atom_types, ' atom types'
+
+        write(11,'(a)') ' '
+        !Write box bd
+        write(11, '(2f23.15, a)') box_bd(1:2), ' xlo xhi'
+        write(11, '(2f23.15, a)') box_bd(3:4), ' ylo yhi'
+        write(11, '(2f23.15, a)') box_bd(5:6), ' zlo zhi'
+
+        !Masses
+        write(11, '(a)') 'Masses'
+
+        write(11, '(a)') ' '
+        do i =1, atom_types
+            call atommass(type_to_name(i),mass)
+            write(11, '(i16, f23.15, 2a)') i, mass, ' # ', type_to_name(i)
+        end do
+        write(11, '(a)') ' '
+
+        write(11, '(a)') 'CAC Elements'
+        write(11, '(a)') ' '
+
+        !Write element nodal positions
+        do i = 1, ele_num
+            select case(trim(adjustl(type_ele(i))))
+            case('fcc')
+                !The first entry is the element specifier 
+                write(11,1) i, basisnum(lat_ele(i)), size_ele(i), size_ele(i), size_ele(i)
+                do ibasis = 1, basisnum(lat_ele(i))
+                    do inod = 1, 8
+                        !Nodal information for every node
+                        write(11,3) inod, ibasis, basis_type(ibasis,lat_ele(i)), r_node(:,ibasis,inod,i)
+                    end do
+                end do
+            end select
+        end do
+
+        do i = 1, atom_num
+            !Element specifier dictating that it is an atom
+            write(11,2) ele_num+i, 1, 1, 1, 1
+            !Write the atomic information
+            write(11,3) 1, 1, type_atom(i), r_atom(:,i)
+        end do
+    end subroutine write_lmpcac
+
     subroutine write_vtk(file)
         !This subroutine writes out a vtk style dump file
         integer :: i, j, inod, ibasis
@@ -257,7 +363,7 @@ module io
         do i = 1, ele_num
             do inod=1, ng_node(lat_ele(i))  
                 do ibasis = 1, basisnum(lat_ele(i))
-                    write(11, '(3f23.1)') sum(r_node(:,:,inod,i),2)/basisnum(lat_ele(i))
+                    write(11, '(3f23.15)') sum(r_node(:,:,inod,i),2)/basisnum(lat_ele(i))
                 end do 
             end do 
         end do
@@ -276,4 +382,29 @@ module io
         end do
         close(11)
     end subroutine
-end module io
\ No newline at end of file
+
+
+! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! READ SUBROUTINES !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!     subroutine read_lmpcac(file, box_bd)
+!         !This subroutine reads in a lmpcac file which can be used with different options and modes
+
+!         !Arguments
+!         character(len=100), intent(in) :: file
+!         real(kind=wp), dimension(6), intent(out) :: box_bd
+
+!         !Internal variables
+!         character(len=1000) :: line
+!         integer :: read_num, atom_lim, ele_lim
+
+!         !Open the lmpcac file
+!         open(unit=11, file=file, action='read', position='rewind')
+
+!         !Skip header lines
+!         read(11,*) line
+!         read(11,*) line
+
+!         !Read total number of elements
+        
+
+!     end subroutine read_lmpcac
+end module io
diff --git a/src/mode_create.f90 b/src/mode_create.f90
index 52342e2..e822798 100644
--- a/src/mode_create.f90
+++ b/src/mode_create.f90
@@ -432,7 +432,6 @@ module mode_create
             end do
 
             !Now figure out how many lattice points could not be contained in elements
-            print *, count(lat_points)
             allocate(r_atom_lat(3,count(lat_points)))
             lat_atom_num = 0
             do ix  = 1, bd_in_array(3)
@@ -453,7 +452,6 @@ module mode_create
                 end do
             end do
 
-            print *, lat_atom_num
         end if 
 
     end subroutine build_with_rhomb

From c1150ae4b13b1418d9419452b455d41d028310fb Mon Sep 17 00:00:00 2001
From: Alex Selimov <alsel@tuta.io>
Date: Fri, 29 Nov 2019 22:23:22 -0500
Subject: [PATCH 02/17] First attempts at figuring out the element format for
 lammpscac

---
 src/io.f90          | 14 ++++++++++----
 src/mode_create.f90 | 46 ++++++++++++++++++++++++++++++++++++++-------
 2 files changed, 49 insertions(+), 11 deletions(-)

diff --git a/src/io.f90 b/src/io.f90
index 948fd52..af66da3 100644
--- a/src/io.f90
+++ b/src/io.f90
@@ -8,6 +8,7 @@ module io
 
     integer :: outfilenum = 0, infilenum = 0
     character(len=100) :: outfiles(10), infiles(10)
+    logical lmpcac
 
    public 
    contains
@@ -57,7 +58,12 @@ module io
                 cycle
             end if
             select case(temp_outfile(scan(temp_outfile,'.',.true.)+1:))
-            case('xyz','lmp','vtk','cac')
+            case('xyz','lmp','vtk')
+                outfilenum=outfilenum+1
+                outfiles(outfilenum) =  temp_outfile
+                exit
+            case('cac')
+                lmpcac = .true. 
                 outfilenum=outfilenum+1
                 outfiles(outfilenum) =  temp_outfile
                 exit
@@ -98,18 +104,18 @@ module io
                 cycle
             end if
 
-            if (scan(temp_outfile,'.',.true.) == 0) then 
+            if (scan(temp_infile,'.',.true.) == 0) then 
                 print *, "No extension included on filename, please type a full filename that includes an extension."
                 read(*,*) temp_infile
                 cycle
             end if
-            select case(temp_outfile(scan(temp_outfile,'.',.true.)+1:))
+            select case(temp_infile(scan(temp_infile,'.',.true.)+1:))
             case('cac')
                 infilenum=infilenum+1
                 infiles(infilenum) =  temp_infile
                 exit
             case default    
-                print *, "File type: ", trim(temp_infile(scan(temp_outfile,'.',.true.):)), " not currently accepted. ", &
+                print *, "File type: ", trim(temp_infile(scan(temp_infile,'.',.true.):)), " not currently accepted. ", &
                          "please input a filename with extension from following list: cac."
                 read(*,*) temp_infile
 
diff --git a/src/mode_create.f90 b/src/mode_create.f90
index e822798..53880c9 100644
--- a/src/mode_create.f90
+++ b/src/mode_create.f90
@@ -11,7 +11,7 @@ module mode_create
 
     character(len=100) :: name, element_type
     real(kind = dp) :: lattice_parameter, orient(3,3), cell_mat(3,8), box_len(3), basis(3,3), origin(3), maxlen(3), &
-                       orient_inv(3,3), box_vert(3,8), maxbd(3), lattice_space(3)
+                       orient_inv(3,3), box_vert(3,8), maxbd(3), lattice_space(3), adjustVar(3,8)
     integer :: esize, duplicate(3), ix, iy, iz, box_lat_vert(3,8), lat_ele_num, lat_atom_num, bd_in_lat(6)
     logical :: dup_flag, dim_flag
 
@@ -91,14 +91,30 @@ module mode_create
             continue
         else
 
-            call cell_init(lattice_parameter, esize, element_type, orient, cell_mat) 
+            if(lmpcac) then 
+                 do inod = 1, 8
+                    do i = 1,3
+                        if(is_equal(cubic_cell(i, inod),0.0_dp)) then 
+                            adjustVar(i,inod) = -0.5_dp*esize
+                        else 
+                            adjustVar(i, inod) = 0.5_dp*esize
+                        end if
+                    end do
+                end do
+            else
+                adjustVar(:,:)=0.0_dp
+            end if
+
+            ! call cell_init(lattice_parameter, esize, element_type, orient, cell_mat) 
             !If the user doesn't pass any build instructions than we just put the cell mat into the element_array 
             call alloc_ele_arrays(1,0)
 
             !Add the basis atoms to the unit cell
             do inod = 1, max_ng_node
                 do ibasis = 1, basisnum(1)
-                    r_node_temp(:,ibasis,inod) = cell_mat(:,inod) + basis_pos(:,ibasis,1) + origin(:)
+                    r_node_temp(:,ibasis,inod) = lattice_parameter*matmul(orient, &
+                                                 matmul(fcc_mat, esize*(cubic_cell(:,inod)+adjustVar(:,inod)))) &
+                                                        + basis_pos(:,ibasis,1)
                 end do
             end do
             do i = 1,3
@@ -276,7 +292,7 @@ module mode_create
         !Internal variables
         integer :: i,  inod, bd_in_lat(6), bd_in_array(6), ix, iy, iz,  numlatpoints, templatpoints, ele(3,8), rzero(3), ilat, &
                    type_interp(basisnum(1)*esize**3), vlat(3), temp_lat(3,8), m, n, o
-        real(kind=dp) :: v(3), temp_nodes(3,1,8), ele_atoms(3,esize**3), r_interp(3,basisnum(1)*esize**3)
+        real(kind=dp) :: v(3), temp_nodes(3,1,8), ele_atoms(3,esize**3), r_interp(3,basisnum(1)*esize**3), adjustVar(3,8)
         real(kind=dp), allocatable :: resize_lat_array(:,:)
         logical, allocatable :: lat_points(:,:,:)
         logical :: node_in_bd(8)
@@ -284,6 +300,21 @@ module mode_create
         !Do some value initialization
         max_esize = esize
 
+        !If we are writing out to lammpscac format then we have to adjust the nodal positions
+
+        if(lmpcac) then 
+             do inod = 1, 8
+                do i = 1,3
+                    if(is_equal(cubic_cell(i, inod),0.0_dp)) then 
+                        adjustVar(i,inod) = -0.25_dp
+                    else 
+                        adjustVar(i, inod) = 0.25_dp
+                    end if
+                end do
+            end do
+        else
+            adjustVar(:,:)=0.0_dp
+        end if
         !First find the bounding lattice points (min and max points for the box in each dimension)
         numlatpoints = 1
         do i = 1, 3
@@ -301,7 +332,6 @@ module mode_create
             continue
         end select
 
-
         !Loop over all of lattice points within the boundary, we choose between two loops. One for the atomistic case 
         !and one for the regular case
         if (esize==2) then  
@@ -415,8 +445,10 @@ module mode_create
 
                         if(all(node_in_bd)) then 
                             lat_ele_num = lat_ele_num+1
-                            r_lat(:,:,lat_ele_num) = temp_nodes(:,1,:)
-
+                            do inod = 1, 8
+                                r_lat(:,inod,lat_ele_num) = temp_nodes(:,1,inod) + adjustVar(:,inod)
+                            end do
+                            
                             !Now set all the lattice points contained within an element to false
                             do o = minval(temp_lat(3,:)), maxval(temp_lat(3,:))
                                 do n = minval(temp_lat(2,:)), maxval(temp_lat(2,:))

From 591762d2591ee26b79313f1ea661685b0b7edf4c Mon Sep 17 00:00:00 2001
From: Alex Selimov <alsel@tuta.io>
Date: Wed, 4 Dec 2019 21:03:13 -0500
Subject: [PATCH 03/17] Latest attempts to get the lammpscac data format
 working

---
 src/mode_create.f90 | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/mode_create.f90 b/src/mode_create.f90
index 53880c9..6c465a3 100644
--- a/src/mode_create.f90
+++ b/src/mode_create.f90
@@ -95,9 +95,9 @@ module mode_create
                  do inod = 1, 8
                     do i = 1,3
                         if(is_equal(cubic_cell(i, inod),0.0_dp)) then 
-                            adjustVar(i,inod) = -0.5_dp*esize
+                            adjustVar(i,inod) = -0.5_dp
                         else 
-                            adjustVar(i, inod) = 0.5_dp*esize
+                            adjustVar(i, inod) = 0.5_dp
                         end if
                     end do
                 end do
@@ -113,7 +113,7 @@ module mode_create
             do inod = 1, max_ng_node
                 do ibasis = 1, basisnum(1)
                     r_node_temp(:,ibasis,inod) = lattice_parameter*matmul(orient, &
-                                                 matmul(fcc_mat, esize*(cubic_cell(:,inod)+adjustVar(:,inod)))) &
+                                                 matmul(fcc_mat, (esize+1)*cubic_cell(:,inod)+adjustVar(:,inod))) &
                                                         + basis_pos(:,ibasis,1)
                 end do
             end do

From b0941e44820615d26889e6eca87ac6b1b5fef55c Mon Sep 17 00:00:00 2001
From: Alex Selimov <alsel@tuta.io>
Date: Tue, 10 Dec 2019 19:33:46 -0500
Subject: [PATCH 04/17] Fixes to file reading to ensure that mode_merge works
 correctly

---
 src/Makefile        |  4 +--
 src/io.f90          | 28 ++++++++++++++-----
 src/mode_create.f90 | 66 ++++++++++++++++++++++++++++++---------------
 3 files changed, 68 insertions(+), 30 deletions(-)

diff --git a/src/Makefile b/src/Makefile
index 22c84eb..c84186e 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
-#FFLAGS=-c -mcmodel=large -Ofast
+#FFLAGS=-mcmodel=large -g -O0 -stand f08 -fpe0 -traceback -check bounds,uninit -warn all -implicitnone -no-wrap-margin
+FFLAGS=-mcmodel=large -Ofast
 MODES=mode_create.o mode_merge.o mode_convert.o
 OBJECTS=main.o elements.o io.o subroutines.o functions.o atoms.o call_mode.o box.o $(MODES)
 
diff --git a/src/io.f90 b/src/io.f90
index 20fac52..e6c4645 100644
--- a/src/io.f90
+++ b/src/io.f90
@@ -394,10 +394,12 @@ module io
         real(kind=dp), dimension(3), intent(in) :: displace
         real(kind = dp), dimension(6), intent(out) :: temp_box_bd
 
-        integer :: i, j, k, n, inod, ibasis, type, size, in_atoms, in_eles
+        integer :: i, j, k, n, inod, ibasis, type, size, in_atoms, in_eles, new_atom_types, &
+                   new_type_to_type(10), new_lattice_types
         character(len=100) :: etype
         real(kind=dp) ::  r(3), newdisplace(3)
         real(kind=dp), allocatable :: r_innode(:,:,:)
+        character(len = 2) :: new_type_to_name(10)
         !First open the file
         open(unit=11, file=trim(adjustl(file)), action='read',position='rewind')
 
@@ -433,15 +435,26 @@ module io
             sub_box_num = sub_box_num + n
 
         !Read in the number of atom types and all their names
-        read(11, *) atom_types, (type_to_name(i), i = 1, atom_types)
+        read(11, *) new_atom_types, (new_type_to_name(i), i = 1, new_atom_types)
+        !Now fit these into the global list of atom types, after this new_type_to_type is the actual global 
+        !type of the atoms within this file
+        do i = 1, new_atom_types
+            call add_atom_type(new_type_to_name(i), new_type_to_type(i))
+        end do
         !Read the number of lattice types, basisnum, and number of nodes for each lattice type
-        read(11,*) lattice_types, (basisnum(i), i = 1, lattice_types), (ng_node(i), i = 1, lattice_types)
+        read(11,*) new_lattice_types, (basisnum(i), i = lattice_types+1, lattice_types+new_lattice_types), &
+                   (ng_node(i), i = lattice_types+1, lattice_types+new_lattice_types)
         !Define max_ng_node and max_basis_num
         max_basisnum = maxval(basisnum)
         max_ng_node = maxval(ng_node)
         !Read the basis atom types for every lattice
-        read(11,*) ((basis_type(i,j), i = 1, basisnum(j)), j = 1, lattice_types)
-
+        read(11,*) ((basis_type(i,j), i = 1, basisnum(j)), j = lattice_types+1, lattice_types+new_lattice_types)
+        !Convert the basis_atom types
+        do j  = lattice_types+1, lattice_types+new_lattice_types
+            do i = 1, basisnum(j)
+                basis_type(i,j) = new_type_to_type(basis_type(i,j))
+            end do 
+        end do 
         !Read number of elements and atoms and allocate arrays
         read(11, *) in_atoms, in_eles
         call grow_ele_arrays(in_eles, in_atoms)
@@ -450,7 +463,7 @@ module io
         !Read the atoms 
         do i = 1, in_atoms
             read(11,*) j, type, r(:)
-            call add_atom(type, r+newdisplace)
+            call add_atom(new_type_to_type(type), r+newdisplace)
         end do
 
         !Read the elements
@@ -462,11 +475,12 @@ module io
                     r_innode(:,ibasis,inod) = r_innode(:, ibasis, inod) + newdisplace
                 end do
             end do
-
+            type = type + lattice_types
             call add_element(etype, size, type, r_innode)
         end do
 
         !Close the file being read
         close(11)
+        lattice_types = lattice_types + new_lattice_types
     end subroutine read_mb
 end module io
\ No newline at end of file
diff --git a/src/mode_create.f90 b/src/mode_create.f90
index b82cc82..08e201a 100644
--- a/src/mode_create.f90
+++ b/src/mode_create.f90
@@ -12,8 +12,8 @@ module mode_create
 
     character(len=100) :: name, element_type
     real(kind = dp) :: lattice_parameter, orient(3,3), cell_mat(3,8), box_len(3), basis(3,3), origin(3), maxlen(3), &
-                       orient_inv(3,3), box_vert(3,8), maxbd(3), lattice_space(3)
-    integer :: esize, duplicate(3), ix, iy, iz, box_lat_vert(3,8), lat_ele_num, lat_atom_num, bd_in_lat(6), &
+                       orient_inv(3,3), box_vert(3,8), maxbd(3), lattice_space(3), duplicate(3)
+    integer :: esize, ix, iy, iz, box_lat_vert(3,8), lat_ele_num, lat_atom_num, bd_in_lat(6), &
                basis_pos(3,10)
     logical :: dup_flag, dim_flag
 
@@ -58,16 +58,15 @@ module mode_create
 
         allocate(r_node_temp(3,max_basisnum,max_ng_node))
 
-        if(dup_flag) then 
+        !Get the inverse orientation matrix we will need later 
+        call matrix_inverse(orient,3,orient_inv)
 
-            !We initialize the cell with a lattice_parameter of 1 because we will add the lattice parameter later
-            call cell_init(1.0_dp, esize, element_type, orient, cell_mat) 
+        if(dup_flag) then 
             
-
+            !Define box vertices
             do i = 1, 8
-                box_vert(:,i) = duplicate(:)*esize*lattice_space(:)*cubic_cell(:,i) + origin(:) 
+                box_vert(:,i) = duplicate(:)*esize*lattice_space(:)*cubic_cell(:,i) + (origin(:)/lattice_parameter) 
             end do
-            call matrix_inverse(orient,3,orient_inv)
             !Now get the rotated box vertex positions in lattice space. Should be integer units
             box_lat_vert = int(matmul(fcc_inv, matmul(orient_inv, box_vert)))+1
             !Find the new maxlen
@@ -76,21 +75,25 @@ module mode_create
                 box_bd(2*i) = maxval(box_vert(i,:)) - 0.25_dp*lattice_space(i)
                 box_bd(2*i-1) = origin(i)-0.25_dp*lattice_space(i)
             end do
-            !and then call the build function with the correct transformation matrix
-            select case(trim(adjustl(element_type)))
-            case('fcc')
 
-                call build_with_rhomb(box_lat_vert, fcc_mat)
-            case default
-                print *, "Element type ", trim(adjustl(element_type)), " not accepted in mode create, please specify a supported ", &
-                         "element type"
-                stop 3
-            end select
 
-            !Now that it is multiply by the lattice parameter
-            box_bd = box_bd*lattice_parameter
+
         else if(dim_flag) then 
-            continue
+            !As a note everything is defined so that the lattice parameter is multiplied in at the end 
+            !so we have to divide all the real Angstroms units by the lattice parameter
+
+            !Define box_vertices
+            do i = 1, 8
+                box_vert(:,i) = (cubic_cell(:,i)*box_len(:) + origin(:))/lattice_parameter
+            end do
+            !Now get the rotated box vertex positions in lattice space. Should be integer units
+            box_lat_vert = int(matmul(fcc_inv, matmul(orient_inv, box_vert)))+1
+
+            !Now get the box_bd in lattice units
+            do i = 1, 3
+                box_bd(2*i) = (box_len(i)+origin(i))/lattice_parameter
+                box_bd(2*i-1) = origin(i)/lattice_parameter
+            end do 
         else
 
             call cell_init(lattice_parameter, esize, element_type, orient, cell_mat) 
@@ -112,6 +115,19 @@ module mode_create
         
         !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
         if(dup_flag.or.dim_flag) then 
+            !Call the build function with the correct transformation matrix
+            select case(trim(adjustl(element_type)))
+            case('fcc')
+
+                call build_with_rhomb(box_lat_vert, fcc_mat)
+            case default
+                print *, "Element type ", trim(adjustl(element_type)), " not accepted in mode create, please specify a supported ", &
+                         "element type"
+                stop 3
+            end select
+            !Now that it is built multiply by the lattice parameter
+            box_bd = box_bd*lattice_parameter
+
             !Allocate variables 
             call alloc_ele_arrays(lat_ele_num, lat_atom_num*basisnum(1))
             if(lat_atom_num > 0) then 
@@ -208,13 +224,21 @@ module mode_create
 
             !If the duplicate command is passed then we extract the information on the new bounds.
             case('duplicate')
+                if(dim_flag) STOP "Both duplicate and dim options cannot be used in mode_create"
                 dup_flag = .true.
                 do i = 1, 3
                     call get_command_argument(arg_pos, textholder)
                     read(textholder, *) duplicate(i)
                     arg_pos = arg_pos + 1
                 end do
-
+            case('dim')
+                if(dup_flag) STOP "Both duplicate and dim options cannot be used in mode_create"
+                dim_flag = .true.
+                do i = 1, 3
+                    call get_command_argument(arg_pos, textholder)
+                    read(textholder, *) box_len(i)
+                    arg_pos = arg_pos + 1
+                end do
             case('origin')
                 do i = 1, 3
                     call get_command_argument(arg_pos, textholder)

From c291ec65b4ce50ee36c0e439bcd0474669a47e78 Mon Sep 17 00:00:00 2001
From: Alex Selimov <alsel@tuta.io>
Date: Fri, 13 Dec 2019 13:44:17 -0500
Subject: [PATCH 05/17] Quick bug fix

---
 src/elements.f90 | 7 +++++--
 1 file changed, 5 insertions(+), 2 deletions(-)

diff --git a/src/elements.f90 b/src/elements.f90
index 28e0e5e..1195e58 100644
--- a/src/elements.f90
+++ b/src/elements.f90
@@ -225,8 +225,11 @@ module elements
         
         exists = .false.
         do i=1, 10
-            if(type == type_to_name(i)) exists = .true.
-            inttype = i
+            if(type == type_to_name(i)) then 
+                exists = .true.
+                inttype = i
+                exit
+            end if
         end do 
 
         if (exists.eqv..false.) then 

From b3a9577cc2e9905cbed879ec6a60316f8628c6b6 Mon Sep 17 00:00:00 2001
From: Alex Selimov <aselimov3@gatech.edu>
Date: Sun, 15 Dec 2019 17:55:31 -0500
Subject: [PATCH 06/17] Added writing to pycac restart file format

---
 src/box.f90      |   9 +++-
 src/elements.f90 |   7 ++-
 src/io.f90       | 136 +++++++++++++++++++++++++++++++++++++++++++----
 3 files changed, 140 insertions(+), 12 deletions(-)

diff --git a/src/box.f90 b/src/box.f90
index ad5fe55..3043851 100644
--- a/src/box.f90
+++ b/src/box.f90
@@ -5,7 +5,7 @@ module box
     implicit none 
 
     real(kind=dp) :: box_bd(6) !Global box boundaries
-
+    character(len=3) :: box_bc !Box boundary conditions (periodic or shrinkwrapped)
     !The subbox variables contain values for each subbox, being the boxes read in through some
     !command. Currently only mode_merge will require sub_boxes, for mode_create it will always 
     !allocate to only 1 sub_box
@@ -14,12 +14,19 @@ module box
     real(kind=dp), allocatable :: sub_box_ori(:,:,:)!Orientations for each of the subboxes
     real(kind=dp), allocatable :: sub_box_bd(:,:)!Boundaries for each of the sub_boxes
 
+
+    !Below are some simulation parameters which may be adjusted if reading in restart files
+    integer :: timestep=0
+    real(kind=dp) :: total_time=0.0_dp
+
+
     public 
     contains 
 
     subroutine box_init
         !Initialize some box functions
         box_bd(:) = 0.0_dp
+        box_bc = 'ppp'
     end subroutine box_init
 
     subroutine alloc_sub_box(n)
diff --git a/src/elements.f90 b/src/elements.f90
index 1195e58..22e6eff 100644
--- a/src/elements.f90
+++ b/src/elements.f90
@@ -11,8 +11,8 @@ module elements
     integer, allocatable :: size_ele(:), lat_ele(:) !Element siz
     real(kind=dp), allocatable :: r_node(:,:,:,:) !Nodal position array
 
-    integer :: ele_num=0 !Number of elements
-    integer :: node_num=0 !Total number of nodes
+    integer, save :: ele_num !Number of elements
+    integer, save :: node_num !Total number of nodes
 
     !Data structure used to represent atoms 
     integer, allocatable ::  type_atom(:)!atom type 
@@ -85,6 +85,9 @@ module elements
         max_basisnum = 0
         basisnum(:) = 0
         ng_node(:) = 0
+        node_num = 0
+        ele_num = 0
+        atom_num = 0
     end subroutine lattice_init
 
     subroutine cell_init(lapa,esize,ele_type, orient_mat, cell_mat)
diff --git a/src/io.f90 b/src/io.f90
index e6c4645..6e914b2 100644
--- a/src/io.f90
+++ b/src/io.f90
@@ -58,7 +58,7 @@ module io
                 cycle
             end if
             select case(temp_outfile(scan(temp_outfile,'.',.true.)+1:))
-            case('xyz', 'lmp', 'vtk', 'mb')
+            case('xyz', 'lmp', 'vtk', 'mb', 'restart')
                 outfilenum=outfilenum+1
                 outfiles(outfilenum) =  temp_outfile
                 exit
@@ -92,6 +92,8 @@ module io
                 call write_vtk(outfiles(i))
             case('mb')
                 call write_mb(outfiles(i))
+            case('restart')
+                call write_pycac(outfiles(i))
             case default
                 print *, "The extension ", trim(adjustl(outfiles(i)(scan(outfiles(i),'.',.true.)+1:))), &
                          " is not accepted for writing. Please select from: xyz and try again"
@@ -107,16 +109,10 @@ module io
         !This is the simplest visualization subroutine, it writes out all nodal positions and atom positions to an xyz file
         character(len=100), intent(in) :: file 
 
-        integer :: node_num, i, inod, ibasis
+        integer :: i, inod, ibasis
 
         open(unit=11, file=trim(adjustl(file)), action='write', status='replace',position='rewind')
 
-        !Calculate total node number
-        node_num=0
-        do i = 1, ele_num
-            node_num = node_num + basisnum(lat_ele(i))*ng_node(lat_ele(i))
-        end do
-
         !Write total number of atoms + elements
         write(11, '(i16)') node_num+atom_num
          
@@ -273,6 +269,123 @@ module io
         close(11)
     end subroutine
 
+    subroutine write_pycac(file)
+        !This subroutine writes restart files meant to be used with the McDowell Group CAC code.
+        !NOTE: This code doesn't work for arbitrary number of basis atoms per node. It assumes that the 
+        !each element has only 1 atom type at the node.
+        character(len=100), intent(in) :: file
+        integer :: interp_max, i, j,  lat_size, inod, ibasis, ip
+        real(kind=dp) :: box_vec(3)
+
+1  format('time' / i16, f23.15)
+2  format('number of elements' / i16)
+3  format('number of nodes' / i16)
+4  format('element types' / i16)
+5  format('number of atoms' / i16)
+6  format('number of grains' / i16)
+7  format('boundary ' / 3a1)
+8  format('box bound' / 6f23.15)
+9  format('box length' / 3f23.15)
+10 format('box matrix')
+11 format(3f23.15)
+12 format('coarse-grained domain')
+13 format('ie ele_type grain_ele lat_type_ele'/ 'ip ibasis type x y z')
+14 format('atomistic domain' / 'ia grain_atom type_atom x y z')
+15 format('maximum lattice periodicity length' / 3f23.15)
+16 format('Number of lattice types and atom types '/ 2i16)
+17 format('lattice type IDs')
+18 format('lattice types for grains')
+19 format('max nodes per element' / i16)
+20 format('max interpo per element' / i16)
+21 format('atom types to elements')
+
+        open(unit=11, file=trim(adjustl(file)), action='write', status='replace',position='rewind')
+
+        write(11,1) timestep, total_time
+        write(11,2) ele_num
+
+        !Below writes the header information for the restart file
+
+        !Calculate the max number of atoms per element
+        select case(max_ng_node)
+        case(8)
+            interp_max = (max_esize)**3
+        end select
+        write(11,20) interp_max
+        write(11,3) node_num
+        write(11,19) max_ng_node
+        write(11,4) lattice_types
+        write(11,2) atom_num
+        write(11,6) 1 !Grain_num is ignored 
+        write(11,16) lattice_types, atom_types
+        write(11,21)
+        do i = 1, atom_types
+            write(11,*) i, type_to_name(i)
+        end do
+        write(11,7) box_bc(1:1), box_bc(2:2), box_bc(3:3)
+        write(11,18)
+        write(11,'(2i16)') 1,1 !This is another throwaway line that is meaningless
+        write(11,17)
+        !This may have to be updated in the future but currently the only 8 node element is fcc
+        do i = 1, lattice_types
+            select case(ng_node(i))
+            case(8)
+                write(11, *) i, 'fcc'
+            end select
+        end do
+        write(11,15) 1.0_dp, 1.0_dp, 1.0_dp !Another throwaway line that isn't needed
+        write(11,8) box_bd
+        write(11,9) box_bd(2)-box_bd(1), box_bd(4) - box_bd(3), box_bd(6)-box_bd(5)
+        write(11,10)
+        !Current boxes are limited to being rectangular 
+        do i = 1,3 
+            box_vec(:) = 0.0_dp
+            box_vec(i) = box_bd(2*i) - box_bd(2*i-1)
+            write(11,11) box_vec
+        end do
+        !We write this as box_mat ori and box_mat current
+        do i = 1,3 
+            box_vec(:) = 0.0_dp
+            box_vec(i) = box_bd(2*i) - box_bd(2*i-1)
+            write(11,11) box_vec
+        end do
+
+        !write the element information
+        if(ele_num > 0) then
+            write(11,12)
+            do i = 1, lattice_types
+                do j = 1, ele_num
+                    if (lat_ele(j) == i) then 
+                        lat_size = size_ele(j)-1
+                        exit
+                    end if
+                end do 
+                write(11,'(3i16)') i, lat_size, basis_type(1,i)
+            end do
+            ip = 0
+            write(11,13)
+            do i = 1, ele_num
+                write(11, '(4i16)') i, lat_ele(i), 1, basis_type(1,lat_ele(i))
+                do inod = 1, ng_node(lat_ele(i))
+                    do ibasis = 1, basisnum(lat_ele(i))
+                        ip = ip + 1
+                        write(11, '(2i16, 3f23.15)') ip, ibasis, r_node(:, ibasis, inod, i)
+                    end do
+                end do
+            end do 
+        end if
+
+        !Now write the atomic information
+        if(atom_num /= 0) then 
+            write(11,14)
+            do i = 1, atom_num
+                write(11, '(3i16, 3f23.15)') i, 1, type_atom(i), r_atom(:,i)
+            end do
+        end if
+
+        close(11)
+    end subroutine write_pycac
+
     subroutine write_mb(file)
 
         !This subroutine writes the cacmb formatted file which provides necessary information for building models
@@ -321,8 +434,10 @@ module io
             end do
         end do
 
+        close(11)
     end subroutine write_mb
 
+
     !!!!!!!!!!!!! Below are subroutines for reading files !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
     subroutine get_in_file(filename)
@@ -481,6 +596,9 @@ module io
 
         !Close the file being read
         close(11)
-        lattice_types = lattice_types + new_lattice_types
+
+        !Only increment the lattice types if there are elements, if there are no elements then we 
+        !just overwrite the arrays
+        if(in_eles > 0) lattice_types = lattice_types + new_lattice_types
     end subroutine read_mb
 end module io
\ No newline at end of file

From 0dff3b2ca9d18eb661e5ee99c7b82a0500544c53 Mon Sep 17 00:00:00 2001
From: Alex Selimov <alsel@tuta.io>
Date: Wed, 11 Dec 2019 09:07:59 -0500
Subject: [PATCH 07/17] Restructure the command parsing loops so the mode
 commands only parse the mode options

---
 src/call_mode.f90    | 10 +++++-----
 src/main.f90         | 30 +++++++++++++++++++++++-------
 src/mode_convert.f90 |  8 +++-----
 src/mode_create.f90  | 33 +++++++++------------------------
 src/mode_merge.f90   | 20 +++++++++-----------
 5 files changed, 49 insertions(+), 52 deletions(-)

diff --git a/src/call_mode.f90 b/src/call_mode.f90
index f571520..fb187c1 100644
--- a/src/call_mode.f90
+++ b/src/call_mode.f90
@@ -1,4 +1,4 @@
-subroutine call_mode(arg_num,mode)
+subroutine call_mode(arg_pos,mode)
     !This code is used to parse the command line argument for the mode information and calls the required 
     !mode module.
 
@@ -9,16 +9,16 @@ subroutine call_mode(arg_num,mode)
 
     implicit none    
 
-    integer, intent(in) :: arg_num
+    integer, intent(out) :: arg_pos
     character(len=100), intent(in) :: mode
 
     select case(mode)
     case('--create')
-        call create
+        call create(arg_pos)
     case('--convert')
-        call convert   
+        call convert(arg_pos)  
     case('--merge')
-        call merge 
+        call merge(arg_pos) 
     case default
         print *, "Mode ", trim(adjustl(mode)), " currently not accepted. Please check documentation for ", &
                  "accepted modes and rerun."
diff --git a/src/main.f90 b/src/main.f90
index c7ac008..a681d33 100644
--- a/src/main.f90
+++ b/src/main.f90
@@ -17,8 +17,8 @@ program main
     use io
     
 
-    integer ::  arg_num
-    character(len=100) ::  mode
+    integer ::  i, end_mode_arg, arg_num
+    character(len=100) ::  argument
 
     !Call initialization functions
     call lattice_init
@@ -29,13 +29,29 @@ program main
 
     !Determine if a mode is being used and what it is. The first argument has to be the mode
     !if a mode is being used
-    call get_command_argument(1, mode)
+    call get_command_argument(1, argument)
 
-    mode = trim(adjustl(mode))
-    if (mode(1:2) == '--') then 
-        call call_mode(arg_num, mode)
+    argument = trim(adjustl(argument))
+    if (argument(1:2) == '--') then 
+        call call_mode(end_mode_arg, argument)
     end if
 
-    !Finish by writing the files
+    !Now we loop through all of the arguments and check for passed options or for a filename to be written out
+    do i = end_mode_arg-1, arg_num
+        call get_command_argument(i, argument)
+
+        !Check to see if a filename was passed
+        if(scan(argument,'.',.true.) > 0) then 
+            call get_out_file(argument)
+        end if
+    end do
+
+    !Check to make sure a file was passed to be written out and then write out 
+            ! Before building do a check on the file
+    if (outfilenum == 0) then
+        argument = 'none'
+        call get_out_file(argument)
+    end if
     call write_out
+
 end program main
\ No newline at end of file
diff --git a/src/mode_convert.f90 b/src/mode_convert.f90
index d6e1d9b..1f8a6bd 100644
--- a/src/mode_convert.f90
+++ b/src/mode_convert.f90
@@ -8,8 +8,9 @@ module mode_convert
     public 
     contains
 
-    subroutine convert
+    subroutine convert(arg_pos)
         !This subroutine converts a single input file from one format to another
+        integer, intent(out) :: arg_pos
         character(len=100) :: infile, outfile
         real(kind = dp) :: temp_box_bd(6)
         !We have to allocate the element and atom arrays with a size of 1 for the read in code to work
@@ -19,10 +20,7 @@ module mode_convert
         call get_in_file(infile)
         call read_in(1, (/0.0_dp,0.0_dp,0.0_dp/), temp_box_bd)
         call grow_box(temp_box_bd)
-
-        !Now get the outfile, writing is done after all the codes complete
-        call get_command_argument(3, outfile)
-        call get_out_file(outfile)
+        arg_pos = 3
 
     end subroutine convert
 end module mode_convert
\ No newline at end of file
diff --git a/src/mode_create.f90 b/src/mode_create.f90
index 08e201a..23b9100 100644
--- a/src/mode_create.f90
+++ b/src/mode_create.f90
@@ -21,11 +21,12 @@ module mode_create
     public 
     contains
 
-    subroutine create()
+    subroutine create(arg_pos)
         ! Main subroutine which controls execution
 
         character(len=100) :: textholder
-    
+        integer, intent(out) :: arg_pos
+
         integer :: i, ibasis, inod
         real(kind=dp), allocatable :: r_node_temp(:,:,:)
 
@@ -45,14 +46,8 @@ module mode_create
         lat_atom_num = 0
 
         !First we parse the command
-        call parse_command()
+        call parse_command(arg_pos)
 
-        ! Before building do a check on the file
-        if (outfilenum == 0) then
-            textholder = 'none'
-            call get_out_file(textholder)
-        end if
-        
         !Now we setup the unit element and call other init subroutines
         call def_ng_node(1, element_type)
 
@@ -161,10 +156,10 @@ module mode_create
         sub_box_array_bd(2,2,1) = ele_num
     end subroutine create
     !This subroutine parses the command and pulls out information needed for mode_create
-    subroutine parse_command()
+    subroutine parse_command(arg_pos)
 
-
-        integer :: arg_pos, ori_pos, i, j, arglen, stat
+        integer, intent(out) :: arg_pos
+        integer ::  ori_pos, i, j, arglen, stat
         character(len=100) :: textholder
         character(len=8) :: orient_string
 
@@ -245,20 +240,10 @@ module mode_create
                     read(textholder, *) origin(i)
                     arg_pos = arg_pos + 1
                 end do
-            !If a filetype is passed then we add name.ext to the outfiles list
-            case('xyz')
-                textholder = trim(adjustl(name)) //'.xyz'
-                call get_out_file(textholder)
 
             case default
-                !Check to see if it is an option command, if so then mode_create must be finished
-                if(textholder(1:1) == '-') then 
-                    exit
-
-                !Check to see if a filename was passed
-                elseif(scan(textholder,'.',.true.) > 0) then 
-                    call get_out_file(textholder)
-                end if
+                !If it isn't an option then you have to exit
+                exit
             end select
         end do
 
diff --git a/src/mode_merge.f90 b/src/mode_merge.f90
index 0919b70..ca9d068 100644
--- a/src/mode_merge.f90
+++ b/src/mode_merge.f90
@@ -12,13 +12,14 @@ module mode_merge
 
     public 
     contains
-    subroutine merge
+    subroutine merge(arg_pos)
         
+        integer, intent(out) :: arg_pos 
         integer :: i
         real(kind=dp) :: displace(3), temp_box_bd(6)
 
         !First we parse the merge command
-        call parse_command
+        call parse_command(arg_pos)
 
         !Now loop over all files and stack them 
         do i = 1, in_num
@@ -44,10 +45,12 @@ module mode_merge
         return
     end subroutine merge
 
-    subroutine parse_command
+    subroutine parse_command(arg_pos)
+
+        integer, intent(out) :: arg_pos
 
         character(len=100) :: textholder
-        integer :: i, stat, arglen, arg_pos
+        integer :: i, stat, arglen
 
         !Get dimension to concatenate along
         call get_command_argument(2, dim, arglen)
@@ -88,13 +91,8 @@ module mode_merge
             select case(trim(textholder))
 
             case default
-                !Check to see if it is an option command, if so then mode_create must be finished
-                if(textholder(1:1) == '-') then 
-                    exit
-                !Check to see if a filename was passed
-                elseif(scan(textholder,'.',.true.) > 0) then 
-                    call get_out_file(textholder)
-                end if
+                !If it isn't an available option to mode merge then we just exit
+                exit 
             end select
         end do
 

From 41024b96743df63463630027198d4846054b21c1 Mon Sep 17 00:00:00 2001
From: Alex Selimov <aselimov3@gatech.edu>
Date: Sun, 15 Dec 2019 18:59:21 -0500
Subject: [PATCH 08/17] Added new parse ori_vec subroutine to parse orientation
 vectors

---
 src/mode_create.f90 | 15 +--------------
 src/subroutines.f90 | 25 +++++++++++++++++++++++++
 2 files changed, 26 insertions(+), 14 deletions(-)

diff --git a/src/mode_create.f90 b/src/mode_create.f90
index 23b9100..bac9b12 100644
--- a/src/mode_create.f90
+++ b/src/mode_create.f90
@@ -199,20 +199,7 @@ module mode_create
                     call get_command_argument(arg_pos, orient_string, arglen)
                     if (arglen==0) STOP "Missing orientation in orient command of mode create"
                     arg_pos = arg_pos+1
-                    ori_pos=2
-                    do j = 1,3
-                        if (orient_string(ori_pos:ori_pos) == '-') then 
-                            ori_pos = ori_pos + 1
-                            read(orient_string(ori_pos:ori_pos), *, iostat=stat) orient(i,j)
-                            if (stat>0) STOP "Error reading orient value"
-                            orient(i,j) = -orient(i,j)
-                            ori_pos = ori_pos + 1
-                        else
-                            read(orient_string(ori_pos:ori_pos), *, iostat=stat) orient(i,j)
-                            if(stat>0) STOP "Error reading orient value"
-                            ori_pos=ori_pos + 1
-                        end if
-                    end do
+                    call parse_ori_vec(orient_string, orient(i,:))
                 end do
 
 
diff --git a/src/subroutines.f90 b/src/subroutines.f90
index 0682037..f4447d7 100644
--- a/src/subroutines.f90
+++ b/src/subroutines.f90
@@ -145,4 +145,29 @@ module subroutines
     return
     end subroutine matrix_inverse
 
+    subroutine parse_ori_vec(ori_string, ori_vec)
+        !This subroutine parses a string to vector in the format [ijk]
+        character(len=8), intent(in) :: ori_string
+        real(kind=dp), dimension(3), intent(out) :: ori_vec
+
+        integer :: i, ori_pos, stat
+
+        ori_pos=2
+        do i = 1,3
+            if (ori_string(ori_pos:ori_pos) == '-') then 
+                ori_pos = ori_pos + 1
+                read(ori_string(ori_pos:ori_pos), *, iostat=stat) ori_vec(i)
+                if (stat>0) STOP "Error reading orientation value"
+                ori_vec(i) = -ori_vec(i)
+                ori_pos = ori_pos + 1
+            else
+                read(ori_string(ori_pos:ori_pos), *, iostat=stat) ori_vec(i)
+                if(stat>0) STOP "Error reading orientation value"
+                ori_pos=ori_pos + 1
+            end if
+        end do
+
+        return
+    end subroutine parse_ori_vec
+
 end module subroutines
\ No newline at end of file

From d0e6253d646928397d43fae8752e908d9008dea1 Mon Sep 17 00:00:00 2001
From: Alex <aselimov3@gatech.edu>
Date: Wed, 25 Dec 2019 16:30:18 -0500
Subject: [PATCH 09/17] Working commit for dislocation generation branch

---
 README.md           | 26 +++++++++++++++++++++++++-
 src/Makefile        | 12 +++++++-----
 src/box.f90         | 12 +++++++++---
 src/functions.f90   | 10 ++++++++++
 src/io.f90          | 12 ++++++++----
 src/main.f90        | 19 ++++++++++++++++---
 src/mode_create.f90 |  2 +-
 src/parameters.f90  |  5 +++++
 src/subroutines.f90 | 33 +++++++++++++++++++++++++++++++++
 9 files changed, 114 insertions(+), 17 deletions(-)

diff --git a/README.md b/README.md
index d390177..c8eb22b 100644
--- a/README.md
+++ b/README.md
@@ -111,4 +111,28 @@ This mode merges multiple data files and creates one big simulation cell. The pa
 
 `N` - The number of files which are being read
 
-`dim` - the dimension they are to be stacked along, can be either `x`, `y`, or `z`. If the argument `none` is passed then the cells are just overlaid. Future options will include a delete overlap command.
\ No newline at end of file
+`dim` - the dimension they are to be stacked along, can be either `x`, `y`, or `z`. If the argument `none` is passed then the cells are just overlaid. Future options will include a delete overlap command.
+
+## Options
+
+Options are additional portions of code which have additional functionality. Options are performed in the order that they appear in the argument list and can be added to any mode. If wanting to use strictly options use `--convert` to specify input and output files.
+
+### Option dislgen
+
+```
+-dislgen [ijk] [hkl] x y z char_angle poisson 
+```
+
+
+
+This options adds an arbitrarily oriented dislocation into your model based on user inputs using the volterra displacement fields. The options are below
+
+`[ijk]` - The vector for the line direction
+
+`[hkl]` - The vector for the slip plane 
+
+`x y z` - The position of the dislocation centroid
+
+`char_angle` - Character angle of the dislocation (0 is screw and 90 is edge)
+
+`poisson` - Poisson's ratio used for the displacement field.
\ No newline at end of file
diff --git a/src/Makefile b/src/Makefile
index c84186e..a55e152 100644
--- a/src/Makefile
+++ b/src/Makefile
@@ -1,8 +1,9 @@
 FC=ifort
 #FFLAGS=-mcmodel=large -g -O0 -stand f08 -fpe0 -traceback -check bounds,uninit -warn all -implicitnone -no-wrap-margin
-FFLAGS=-mcmodel=large -Ofast
+FFLAGS=-mcmodel=large -Ofast -no-wrap-margin
 MODES=mode_create.o mode_merge.o mode_convert.o
-OBJECTS=main.o elements.o io.o subroutines.o functions.o atoms.o call_mode.o box.o $(MODES)
+OPTIONS=opt_disl.o
+OBJECTS=main.o elements.o io.o subroutines.o functions.o atoms.o call_mode.o box.o call_option.o $(MODES) $(OPTIONS)
 
 .SUFFIXES:
 .SUFFIXES: .c .f .f90 .F90 .o
@@ -25,9 +26,10 @@ cleantest:
 	$(RM) testfuncs testfuncs.o
 
 $(OBJECTS) : parameters.o
-atoms.o subroutines.o testfuncs.o : functions.o
-main.o io.o build_subroutines.o: elements.o 
+atoms.o subroutines.o testfuncs.o box.o : functions.o
+main.o io.o $(MODES) $(OPTIONS) : elements.o 
 call_mode.o : $(MODES)
+call_option.o : $(OPTIONS)
 $(MODES) io.o: atoms.o box.o
 $(MODES) main.o : io.o
-testfuncs.o elements.o mode_create.o: subroutines.o
+testfuncs.o elements.o mode_create.o opt_disl.o: subroutines.o
diff --git a/src/box.f90 b/src/box.f90
index 3043851..de34e23 100644
--- a/src/box.f90
+++ b/src/box.f90
@@ -1,7 +1,7 @@
 module box
     !This module contains information on the properties of the current box.
     use parameters
-
+    use functions
     implicit none 
 
     real(kind=dp) :: box_bd(6) !Global box boundaries
@@ -34,8 +34,14 @@ module box
 
         integer, intent(in) :: n
 
-        allocate(sub_box_ori(3,3,n), sub_box_bd(6,n), sub_box_array_bd(2,2,n))
+        integer :: i
 
+        allocate(sub_box_ori(3,3,n), sub_box_bd(6,n), sub_box_array_bd(2,2,n))
+        do i = 1, n
+            sub_box_ori(:,:,i) = identity_mat(3)
+            sub_box_bd(:,i) = 0.0_dp
+            sub_box_array_bd(:,:,i) = 1
+        end do
     end subroutine alloc_sub_box
 
     subroutine grow_sub_box(n)
@@ -58,7 +64,7 @@ module box
         call move_alloc(temp_bd, sub_box_bd)
 
         temp_array_bd(:,:,1:sub_box_num) = sub_box_array_bd
-        temp_array_bd(:,:,sub_box_num+1:) = 0.0_dp
+        temp_array_bd(:,:,sub_box_num+1:) = 1
         call move_alloc(temp_array_bd, sub_box_array_bd)
 
         return
diff --git a/src/functions.f90 b/src/functions.f90
index eb8e813..f340b86 100644
--- a/src/functions.f90
+++ b/src/functions.f90
@@ -232,4 +232,14 @@ END FUNCTION StrDnCase
         end if
     return
     end function is_equal
+
+    pure function unitvec(n,vec)
+        integer, intent(in) :: n
+        real(kind=dp), intent(in) :: vec(n)
+        real(kind=dp) :: unitvec(n)
+
+        unitvec = vec/norm2(vec)
+        return
+    end function unitvec
+
 end module functions
diff --git a/src/io.f90 b/src/io.f90
index 6e914b2..46ece37 100644
--- a/src/io.f90
+++ b/src/io.f90
@@ -409,7 +409,7 @@ module io
         end do
 
         !Write the number of atom types in the current model and all of their names
-        write(11,*) atom_types, (type_to_name(i), i=1, atom_types)
+        write(11,*) atom_types, (type_to_name(i)//' ', i=1, atom_types)
         !Write the number of lattice_types, basisnum and number of nodes for each lattice type
         write(11,*) lattice_types, (basisnum(i), i = 1, lattice_types), (ng_node(i), i = 1, lattice_types)
         !Now for every lattice type write the basis atom types
@@ -544,10 +544,14 @@ module io
             read(11,*) sub_box_bd(:,sub_box_num+i)
             sub_box_bd(:,sub_box_num+i) = sub_box_bd(:, sub_box_num+i) + displace(:)
             !Read in sub_box_array_bd
-            read(11,*) ((sub_box_ori(j, k, sub_box_num+i), j = 1, 2), k = 1, 2)
-
+            read(11,*) ((sub_box_array_bd(j, k, sub_box_num+i), j = 1, 2), k = 1, 2)
         end do
-            sub_box_num = sub_box_num + n
+
+        !Add the existing element boundaries
+        sub_box_array_bd(:,1,sub_box_num+1:) = sub_box_array_bd(:,1,sub_box_num+1:) + atom_num
+        sub_box_array_bd(:,2,sub_box_num+1:) = sub_box_array_bd(:,2,sub_box_num+1:) + ele_num
+
+        sub_box_num = sub_box_num + n
 
         !Read in the number of atom types and all their names
         read(11, *) new_atom_types, (new_type_to_name(i), i = 1, new_atom_types)
diff --git a/src/main.f90 b/src/main.f90
index a681d33..6bdc0ed 100644
--- a/src/main.f90
+++ b/src/main.f90
@@ -17,7 +17,7 @@ program main
     use io
     
 
-    integer ::  i, end_mode_arg, arg_num
+    integer ::  i, end_mode_arg, arg_num, arg_pos
     character(len=100) ::  argument
 
     !Call initialization functions
@@ -37,12 +37,25 @@ program main
     end if
 
     !Now we loop through all of the arguments and check for passed options or for a filename to be written out
-    do i = end_mode_arg-1, arg_num
-        call get_command_argument(i, argument)
+    arg_pos = end_mode_arg
+    do while(.true.)
+        !Exit the loop if we are done reading
+        if(arg_pos > arg_num) exit
+
+        call get_command_argument(arg_pos, argument)
 
         !Check to see if a filename was passed
         if(scan(argument,'.',.true.) > 0) then 
             call get_out_file(argument)
+            arg_pos = arg_pos + 1
+
+        !Check to see if an option has been passed
+        else if(argument(1:1) == '-') then 
+            call call_option(argument, arg_pos)
+        !Otherwise print that the argument is not accepted and move on
+        else
+            print *, trim(adjustl(argument)), " is not accepted. Skipping to next argument"
+            arg_pos = arg_pos + 1
         end if
     end do
 
diff --git a/src/mode_create.f90 b/src/mode_create.f90
index bac9b12..ee35488 100644
--- a/src/mode_create.f90
+++ b/src/mode_create.f90
@@ -230,10 +230,10 @@ module mode_create
 
             case default
                 !If it isn't an option then you have to exit
+                arg_pos = arg_pos -1
                 exit
             end select
         end do
-
         !Calculate the lattice periodicity length in lattice units
         do i = 1, 3
             lattice_space(i) = norm2(orient(i,:))
diff --git a/src/parameters.f90 b/src/parameters.f90
index 0443622..b54a702 100644
--- a/src/parameters.f90
+++ b/src/parameters.f90
@@ -2,7 +2,12 @@ module parameters
 
     implicit none
 
+    !Default precision
     integer, parameter :: dp= selected_real_kind(15,307)
+    !Parameters for floating point tolerance
     real(kind=dp), parameter :: lim_zero = epsilon(1.0_dp), &
                                 lim_large = huge(1.0_dp)
+    !Numeric constants
+    real(kind=dp), parameter :: pi = 3.14159265358979323846_dp
+
 end module parameters
\ No newline at end of file
diff --git a/src/subroutines.f90 b/src/subroutines.f90
index f4447d7..d0b004b 100644
--- a/src/subroutines.f90
+++ b/src/subroutines.f90
@@ -1,6 +1,7 @@
 module subroutines
     use parameters
     use functions
+    use box
     implicit none
 
     integer :: allostat, deallostat
@@ -170,4 +171,36 @@ module subroutines
         return
     end subroutine parse_ori_vec
 
+    subroutine parse_pos(i, pos_string, pos)
+        !This subroutine parses the pos command allowing for command which include inf
+        integer, intent(in) ::  i !The dimension of the position
+        character(len=100), intent(in) :: pos_string !The position string
+        real(kind=dp), intent(out) :: pos !The output parsed position value
+
+        integer :: iospara
+        if(trim(adjustl(pos_string)) == 'inf') then 
+            pos=box_bd(2*i)
+        else if(trim(adjustl(pos_string)) == '-inf') then
+            pos=box_bd(2*i-1)
+        else if ((index(pos_string,'-') > 0).and.(index(pos_string,'inf')>0)) then 
+            !Now extract the number we are reducing from infinity
+            read(pos_string(index(pos_string,'-')+1:), *, iostat=iospara) pos
+            pos = box_bd(2*i) - pos
+        else if ((index(pos_string,'+') > 0).and.(index(pos_string,'inf')>0)) then 
+            !Now extract the number we are reducing from infinity
+            read(pos_string(index(pos_string,'+')+1:), *, iostat=iospara) pos
+            pos = box_bd(2*i-1) + pos
+        else if ((index(pos_string,'*') > 0).and.(index(pos_string,'inf')>0)) then 
+            !Now extract the number we are reducing from infinity
+            read(pos_string(index(pos_string,'*')+1:), *, iostat=iospara) pos
+            pos = (box_bd(2*i)-box_bd(2*i-1))*pos
+        else
+            read(pos_string, *, iostat=iospara) pos
+        end if
+
+        if (iospara > 0) then 
+            print *, "Error reading position argument ", trim(adjustl(pos_string)), ". Please reformat and try again."
+        end if
+    end subroutine parse_pos
+
 end module subroutines
\ No newline at end of file

From 8fc91f4dd208f4993409ec27a69f68b1b66a5728 Mon Sep 17 00:00:00 2001
From: Alex Selimov <alsel@tuta.io>
Date: Mon, 13 Jan 2020 09:42:57 -0500
Subject: [PATCH 10/17] Working planar vacancy cluster algorithm

---
 README.md           | 16 +++++++++++++++-
 src/Makefile        | 10 +++++++---
 src/box.f90         | 17 +++++++++++++++++
 src/elements.f90    | 23 +++++++++++++++++++++++
 src/functions.f90   |  8 ++++++++
 src/subroutines.f90 | 46 +++++++++++++++++++++++++++++++++++++++++++++
 6 files changed, 116 insertions(+), 4 deletions(-)

diff --git a/README.md b/README.md
index c8eb22b..b4b690d 100644
--- a/README.md
+++ b/README.md
@@ -135,4 +135,18 @@ This options adds an arbitrarily oriented dislocation into your model based on u
 
 `char_angle` - Character angle of the dislocation (0 is screw and 90 is edge)
 
-`poisson` - Poisson's ratio used for the displacement field.
\ No newline at end of file
+`poisson` - Poisson's ratio used for the displacement field.
+
+### Option disloop
+
+````
+-disloop loop_normal loop_size x y z
+````
+
+This option deletes vacancies on a plane which when minimized should result in a dislocation loop structure. The arguments are below:
+
+`dim` - The box dimension which defines the normal to the loop plane. As of now this dimension must be a closed back direction, meaning that for fcc a box dimension has to be of the (111) family of planes
+
+`n` - The number of atoms to delete on the loop plane
+
+`x y z` - The centroid of the loop. As a note, 
\ No newline at end of file
diff --git a/src/Makefile b/src/Makefile
index a55e152..acc2066 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
-FFLAGS=-mcmodel=large -Ofast -no-wrap-margin
+FFLAGS=-mcmodel=large -g -O0 -stand f08 -fpe0 -traceback -check bounds,uninit -warn all -implicitnone -no-wrap-margin
+#FFLAGS=-mcmodel=large -Ofast -no-wrap-margin
 MODES=mode_create.o mode_merge.o mode_convert.o
 OPTIONS=opt_disl.o
 OBJECTS=main.o elements.o io.o subroutines.o functions.o atoms.o call_mode.o box.o call_option.o $(MODES) $(OPTIONS)
@@ -19,12 +19,16 @@ clean:
 	$(RM) cacmb  *.o 
 
 testfuncs: testfuncs.o functions.o subroutines.o
-	$(FC) testfuncs.o functions.o subroutines.o elements.o -o $@
+	$(FC) testfuncs.o functions.o subroutines.o box.o elements.o -o $@
 
 .PHONY: cleantest
 cleantest:
 	$(RM) testfuncs testfuncs.o
 
+.PHONY: test
+test: testfuncs
+	./testfuncs
+
 $(OBJECTS) : parameters.o
 atoms.o subroutines.o testfuncs.o box.o : functions.o
 main.o io.o $(MODES) $(OPTIONS) : elements.o 
diff --git a/src/box.f90 b/src/box.f90
index de34e23..de8ef59 100644
--- a/src/box.f90
+++ b/src/box.f90
@@ -85,4 +85,21 @@ module box
         return
     end subroutine grow_box
 
+
+    subroutine in_sub_box(r, which_sub_box)
+        !This returns which sub_box a point is in. It returns the first sub_box with boundaries which 
+        !contain the point.
+        real(kind=dp), dimension(3), intent(in) :: r
+        integer, intent(out) :: which_sub_box
+
+        integer :: i 
+
+        do i = 1, sub_box_num
+            if( in_block_bd(r, sub_box_bd(:,i))) then 
+                which_sub_box = i
+                exit
+            end if
+        end do 
+    return 
+    end subroutine in_sub_box
 end module box
\ No newline at end of file
diff --git a/src/elements.f90 b/src/elements.f90
index 22e6eff..1c17643 100644
--- a/src/elements.f90
+++ b/src/elements.f90
@@ -3,6 +3,7 @@ module elements
     !This module contains the elements data structures, structures needed for building regions
     !and operations that are done on elements
     use parameters
+    use functions
     use subroutines
     implicit none
 
@@ -355,4 +356,26 @@ module elements
         return
     end subroutine rhombshape
 
+    subroutine delete_atoms(num, index)
+        !This subroutine deletes atoms from the atom arrays 
+        integer, intent(in) :: num
+        integer, intent(inout), dimension(num) :: index
+
+        integer :: i, j
+
+        call heapsort(index)
+
+        !We go from largest index to smallest index just to make sure that we don't miss 
+        !accidentally overwrite values which need to be deleted
+        do i = num, 1, -1
+            if(index(i) == atom_num) then 
+                r_atom(:,index(i)) = 0.0_dp
+                type_atom(index(i)) = 0
+            else
+                r_atom(:,index(i)) = r_atom(:, atom_num)
+                type_atom(index(i)) = type_atom(atom_num)
+            end if            
+            atom_num = atom_num - 1
+        end do 
+    end subroutine
 end module elements
\ No newline at end of file
diff --git a/src/functions.f90 b/src/functions.f90
index f340b86..82bcbcd 100644
--- a/src/functions.f90
+++ b/src/functions.f90
@@ -242,4 +242,12 @@ END FUNCTION StrDnCase
         return
     end function unitvec
 
+    pure function norm_dis(rl,rk)
+        !This just returns the magnitude of the vector between two points
+        real(kind=dp), dimension(3), intent(in) :: rl, rk
+        real(kind=dp) :: norm_dis(4)
+
+        norm_dis(1:3) = (rk - rl)
+        norm_dis(4) = norm2(rk-rl)
+    end function
 end module functions
diff --git a/src/subroutines.f90 b/src/subroutines.f90
index d0b004b..330fbca 100644
--- a/src/subroutines.f90
+++ b/src/subroutines.f90
@@ -203,4 +203,50 @@ module subroutines
         end if
     end subroutine parse_pos
 
+    subroutine heapsort(a)
+     
+       integer, intent(in out) :: a(0:)
+       integer :: start, n, bottom
+       real :: temp
+     
+       n = size(a)
+       do start = (n - 2) / 2, 0, -1
+         call siftdown(a, start, n);
+       end do
+     
+       do bottom = n - 1, 1, -1
+         temp = a(0)
+         a(0) = a(bottom)
+         a(bottom) = temp;
+         call siftdown(a, 0, bottom)
+       end do
+     
+    end subroutine heapsort
+     
+    subroutine siftdown(a, start, bottom)
+     
+      integer, intent(in out) :: a(0:)
+      integer, intent(in) :: start, bottom
+      integer :: child, root
+      real :: temp
+     
+      root = start
+      do while(root*2 + 1 < bottom)
+        child = root * 2 + 1
+     
+        if (child + 1 < bottom) then
+          if (a(child) < a(child+1)) child = child + 1
+        end if
+     
+        if (a(root) < a(child)) then
+          temp = a(child)
+          a(child) = a (root)
+          a(root) = temp
+          root = child
+        else
+          return
+        end if  
+      end do      
+     
+    end subroutine siftdown
 end module subroutines
\ No newline at end of file

From 1f1598dec9f4e4cb0178592b97b23f53914ee0e2 Mon Sep 17 00:00:00 2001
From: Alex Selimov <alsel@tuta.io>
Date: Mon, 13 Jan 2020 09:49:42 -0500
Subject: [PATCH 11/17] Added Make install and fix to main

---
 src/Makefile | 4 ++++
 src/main.f90 | 5 +++++
 2 files changed, 9 insertions(+)

diff --git a/src/Makefile b/src/Makefile
index acc2066..c20e5b2 100644
--- a/src/Makefile
+++ b/src/Makefile
@@ -29,6 +29,10 @@ cleantest:
 test: testfuncs
 	./testfuncs
 
+.PHONY: install
+install: cacmb
+	cp ./cacmb /usr/local/bin
+
 $(OBJECTS) : parameters.o
 atoms.o subroutines.o testfuncs.o box.o : functions.o
 main.o io.o $(MODES) $(OPTIONS) : elements.o 
diff --git a/src/main.f90 b/src/main.f90
index 6bdc0ed..49b1a36 100644
--- a/src/main.f90
+++ b/src/main.f90
@@ -23,6 +23,7 @@ program main
     !Call initialization functions
     call lattice_init
     call box_init
+    end_mode_arg = 0
     
     ! Command line parsing
     arg_num = command_argument_count()
@@ -36,6 +37,10 @@ program main
         call call_mode(end_mode_arg, argument)
     end if
 
+    !Check to make sure a mode was called
+    if (end_mode_arg=0) then 
+        stop "Please run cacmb using an available mode"
+    end if
     !Now we loop through all of the arguments and check for passed options or for a filename to be written out
     arg_pos = end_mode_arg
     do while(.true.)

From ec7980ff0b2186adf813172f305ba4e5f3e1aa2d Mon Sep 17 00:00:00 2001
From: Alex Selimov <alsel@tuta.io>
Date: Mon, 13 Jan 2020 09:50:29 -0500
Subject: [PATCH 12/17] Fix to main when running without a mode

---
 src/main.f90 | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/main.f90 b/src/main.f90
index 49b1a36..ed34c53 100644
--- a/src/main.f90
+++ b/src/main.f90
@@ -38,8 +38,8 @@ program main
     end if
 
     !Check to make sure a mode was called
-    if (end_mode_arg=0) then 
-        stop "Please run cacmb using an available mode"
+    if (end_mode_arg==0) then 
+        stop "Nothing to do, please run cacmb using an available mode"
     end if
     !Now we loop through all of the arguments and check for passed options or for a filename to be written out
     arg_pos = end_mode_arg

From 20546dc2672435e6c0b8c73efb4bd60de47affe8 Mon Sep 17 00:00:00 2001
From: Alex Selimov <alsel@tuta.io>
Date: Mon, 13 Jan 2020 16:43:49 -0500
Subject: [PATCH 13/17] Actual dislocation loop generation algorithm alongside
 dislgen command

---
 README.md        |  10 +-
 src/io.f90       |   2 +-
 src/opt_disl.f90 | 520 +++++++++++++++++++++++++++++++++++++++++++++++
 3 files changed, 528 insertions(+), 4 deletions(-)
 create mode 100644 src/opt_disl.f90

diff --git a/README.md b/README.md
index b4b690d..cf7fcfb 100644
--- a/README.md
+++ b/README.md
@@ -140,13 +140,17 @@ This options adds an arbitrarily oriented dislocation into your model based on u
 ### Option disloop
 
 ````
--disloop loop_normal loop_size x y z
+-disloop loop_normal radius x y z bx by bz poisson
 ````
 
 This option deletes vacancies on a plane which when minimized should result in a dislocation loop structure. The arguments are below:
 
-`dim` - The box dimension which defines the normal to the loop plane. As of now this dimension must be a closed back direction, meaning that for fcc a box dimension has to be of the (111) family of planes
+`dim` - The box dimension which defines the normal to the loop plane. As of now this dimension must be a closed back direction, meaning that for fcc a box dimension has to be of the (111) family of planes. Either `x`, `y`, or `z`.
 
 `n` - The number of atoms to delete on the loop plane
 
-`x y z` - The centroid of the loop. As a note, 
\ No newline at end of file
+`x y z` - The centroid of the loop. 
+
+`bx by bz` - The burgers vector for the dislocation
+
+`poisson` - Poisson ratio for continuum solution
\ No newline at end of file
diff --git a/src/io.f90 b/src/io.f90
index 46ece37..ba26062 100644
--- a/src/io.f90
+++ b/src/io.f90
@@ -315,7 +315,7 @@ module io
         write(11,3) node_num
         write(11,19) max_ng_node
         write(11,4) lattice_types
-        write(11,2) atom_num
+        write(11,5) atom_num
         write(11,6) 1 !Grain_num is ignored 
         write(11,16) lattice_types, atom_types
         write(11,21)
diff --git a/src/opt_disl.f90 b/src/opt_disl.f90
new file mode 100644
index 0000000..4510ede
--- /dev/null
+++ b/src/opt_disl.f90
@@ -0,0 +1,520 @@
+module opt_disl
+
+    !This module contains all code associated with dislocations
+
+    use parameters
+    use elements
+    use subroutines
+    use box 
+    implicit none
+
+
+    real(kind=dp), dimension(3) :: line, slip_plane, centroid!dislocation line, slip plane vectors, centroid,
+    real(kind=dp) ::  burgers(3)  !burgers vector of loop
+    real(kind=dp) :: poisson, char_angle, lattice_parameter!Poisson ratio and character angle, lattice_parameter for burgers vector
+    character(len=10) :: lattice
+    character(len=1) :: loop_normal
+    real(kind=dp) :: loop_radius !Dislocation loop radius
+    
+    real(kind=dp) :: b !burgers vector
+
+
+    public 
+    contains 
+
+    subroutine dislocation(option, arg_pos)
+
+        !Main calling function for all codes related to dislocations
+
+        character(len=100), intent(in) :: option
+        integer, intent(inout) :: arg_pos 
+
+
+        select case(trim(adjustl(option)))
+        case('-dislgen')
+            call parse_dislgen(arg_pos)
+            call dislgen
+        case('-disloop')
+            call parse_disloop(arg_pos)
+            call disloop
+        end select
+    end subroutine dislocation
+
+    subroutine parse_dislgen(arg_pos)
+
+        !Parse dislgen command
+
+        integer, intent(inout) :: arg_pos
+
+        integer :: i,arglen
+        character(len=8) :: ori_string
+        character(len=100) :: textholder
+
+        !Parse all of the commands
+        arg_pos = arg_pos + 1
+        line(:) = 0.0_dp
+        call get_command_argument(arg_pos, ori_string, arglen)
+        if (arglen==0) STOP "Missing line vector in dislgen command"
+        call parse_ori_vec(ori_string, line)
+
+        arg_pos=arg_pos + 1
+        slip_plane(:) = 0.0_dp
+        call get_command_argument(arg_pos, ori_string, arglen)
+        if (arglen==0) STOP "Missing plane vector in dislgen command"
+        call parse_ori_vec(ori_string, slip_plane)
+
+        do i = 1, 3
+            arg_pos = arg_pos + 1
+            call get_command_argument(arg_pos, textholder, arglen)
+            if (arglen==0) STOP "Missing centroid in dislgen command"
+            call parse_pos(i, textholder, centroid(i))           
+        end do
+
+        print *, centroid
+
+        arg_pos = arg_pos + 1
+        call get_command_argument(arg_pos, textholder, arglen)
+        if (arglen==0) STOP "Missing character angle in dislgen command"
+        read(textholder, *) char_angle
+
+
+        arg_pos = arg_pos + 1
+        call get_command_argument(arg_pos, textholder, arglen)
+        if (arglen==0) STOP "Missing poisson in dislgen command"
+        read(textholder, *) poisson
+
+        arg_pos = arg_pos + 1
+        call get_command_argument(arg_pos, lattice, arglen)
+        if (arglen==0) STOP "Missing lattice in dislgen command"
+
+        arg_pos = arg_pos + 1
+        call get_command_argument(arg_pos, textholder, arglen)
+        if (arglen==0) STOP "Missing lattice parameter in dislgen command"
+        read(textholder, *) lattice_parameter
+
+        arg_pos = arg_pos + 1
+        !Now set the vurgers vector based on the lattice paarameter and lattice type
+        select case(lattice)
+            case('fcc')
+                b = lattice_parameter / sqrt(2.0_dp)
+            ! case('bcc')
+            !     b = lattice_parameter * sqrt(3.0_dp) / 2.0_dp
+            case default
+                print *, 'Error: Lattice structure', lattice, ' is not accepted for dislgen option'
+                STOP
+       end select
+
+    end subroutine parse_dislgen
+
+    subroutine dislgen
+        !This subroutine creates the actual dislocation
+
+        integer :: i, sub_box, inod, ibasis
+        real(kind=dp) :: ss_ori(3,3), ss_inv(3,3), be, bs, slipx(3), disp_transform(3,3), inv_transform(3,3), &
+                         actan, r(3), disp(3)
+
+        !Calculate screw and edge burgers vectors
+        be = sin(char_angle*pi/180.0_dp)*b
+        bs = cos(char_angle*pi/180.0_dp)*b
+
+        !Figure out which sub box you are in so you can access it's orientation, this code will not work
+        !with overlapping sub_boxes
+        do i = 1, sub_box_num
+            if(in_block_bd(centroid, sub_box_bd(:,i))) then 
+                sub_box = i
+                exit
+            end if
+        end do
+
+        !Construct the slip system orientation matrix in an unrotated system 
+        slipx = cross_product(slip_plane, line)
+        ss_ori(1,:) = unitvec(3,slipx)
+        ss_ori(2,:) = unitvec(3,slip_plane)
+        ss_ori(3,:) = unitvec(3,line)
+        call matrix_inverse(ss_ori, 3, ss_inv)
+
+        !Apply the rotation
+        disp_transform = matmul(sub_box_ori(:,:,i), ss_inv)
+        call matrix_inverse(disp_transform,3,inv_transform)
+
+        if(atom_num > 0) then 
+            do i = 1, atom_num
+                r=r_atom(:,i) - centroid
+                r=matmul(inv_transform, r)
+                if (r(1) == 0) then 
+                    actan=pi/2
+                else
+                    actan = datan2(r(2),r(1))
+                end if
+
+                if ((r(1)**2 + r(2)**2) == 0) cycle
+
+                !This is the elastic displacement field for dislocation according to Hirth and Lowe
+                disp(1) = be/(2.0_dp*pi) * (actan + (r(1)*r(2))/(2.0_dp*(1.0_dp-poisson)*(r(1)**2.0_dp + r(2)**2.0_dp)))
+                disp(2) = -be/(2.0_dp*pi)*((1.0_dp-2.0_dp*poisson)/(4.0_dp-4.0_dp*poisson) * &
+                                            log(r(1)**2.0_dp + r(2)**2.0_dp) &
+                                            + (r(1)**2.0_dp - r(2)**2.0_dp)/(4.0_dp*(1.0_dp-poisson)&
+                                            *(r(1)**2.0_dp+r(2)**2.0_dp)))
+                disp(3) = bs/(2.0_dp*pi) * actan
+
+                !This transforms the displacement to the correct orientation
+                disp = matmul(disp_transform, disp)
+
+                r_atom(:,i) = r_atom(:,i) + disp
+            end do
+        end if
+
+        if(ele_num > 0) then 
+            do i = 1, ele_num
+                do inod=1, ng_node(lat_ele(i))
+                    do ibasis = 1, basisnum(lat_ele(i))
+                        r = r_node(:,ibasis,inod,i)
+                        r = matmul(inv_transform, r)
+                        if (r(1) == 0) then 
+                            actan = pi/2
+                        else
+                            actan = datan2(r(2),r(1))
+                        end if
+
+                        if ((r(1)**2 + r(2)**2) == 0) cycle
+                
+                        !This is the elastic displacement field for dislocation according to Hirth and Lowe
+                        disp(1) = be/(2.0_dp*pi) * (actan + (r(1)*r(2))/(2.0_dp*(1.0_dp-poisson)*(r(1)**2.0_dp + r(2)**2.0_dp)))
+                        disp(2) = -be/(2.0_dp*pi)*((1.0_dp-2.0_dp*poisson)/(4.0_dp-4.0_dp*poisson) * &
+                                                    log(r(1)**2.0_dp + r(2)**2.0_dp) &
+                                                    + (r(1)**2.0_dp - r(2)**2.0_dp)/(4.0_dp*(1.0_dp-poisson)&
+                                                    *(r(1)**2.0_dp+r(2)**2.0_dp)))
+                        disp(3) = bs/(2.0_dp*pi) * actan
+
+
+                        disp = matmul(disp_transform, disp)
+                        r_node(:,ibasis,inod,i) = r_node(:,ibasis,inod,i) + disp
+                    end do
+                end do
+            end do
+        end if
+
+    end subroutine dislgen
+
+    subroutine parse_disloop(arg_pos)
+        !This subroutine parses the disloop command
+
+        integer, intent(inout) :: arg_pos
+
+        integer :: i,arglen, sbox
+        character(len=8) :: ori_string
+        character(len=100) :: textholder
+
+        !Parse all of the commands
+        arg_pos = arg_pos + 1
+        loop_normal = ' '
+        call get_command_argument(arg_pos, loop_normal, arglen)
+        if (arglen==0) STOP "Missing loop_normal in disloop command"
+        !Convert the loop_normal to the dimension
+        select case(loop_normal)
+            case('x','X', 'y', 'Y', 'z', 'Z')
+               continue
+            case default
+                print *, "Dimension argument must either be x, y, or z not", loop_normal
+                stop 3
+        end select
+
+        arg_pos = arg_pos + 1
+        loop_radius = 0
+        call get_command_argument(arg_pos, textholder, arglen)
+        if (arglen==0) STOP "Missing loop_size in disloop command"
+        read(textholder, *) loop_radius
+
+        do i = 1, 3
+            arg_pos = arg_pos + 1
+            call get_command_argument(arg_pos, textholder, arglen)
+            if (arglen==0) STOP "Missing centroid in disloop command"
+            call parse_pos(i, textholder, centroid(i))           
+        end do
+
+        burgers(:) = 0.0_dp
+        do i = 1, 3
+            arg_pos = arg_pos + 1
+            call get_command_argument(arg_pos, textholder, arglen)
+            if (arglen==0) STOP "Missing burgers vector in disloop command"
+            read(textholder, *) burgers(i)
+        end do
+
+        arg_pos = arg_pos + 1
+        call get_command_argument(arg_pos, textholder, arglen)
+        if (arglen==0) STOP "Missing poisson ratio in disloop command"
+        read(textholder, *) poisson
+
+        arg_pos = arg_pos + 1
+
+        !Now check to make sure that the dimension selected is actually a 1 1 1 direction.
+        ! call in_sub_box(centroid, sbox)
+        ! if(.not.((abs(sub_box_ori(loop_normal,1,sbox)) == abs(sub_box_ori(loop_normal,2,sbox))).and. &
+        !          (abs(sub_box_ori(loop_normal,2,sbox)) == abs(sub_box_ori(loop_normal,3,sbox))).and. &
+        !          (abs(sub_box_ori(loop_normal,3,sbox)) == abs(sub_box_ori(loop_normal,1,sbox))))) then 
+
+        !     print *, "The selected dimension ", loop_normal, " for sub_box ", sbox, " is ", &
+        !              sub_box_ori(loop_normal,:,sbox), " which is not in the (111) family of planes"
+        !     STOP 3 
+        ! end if
+
+
+    end subroutine
+
+    !Code for the creation of dislocation loops is based on functions from atomsk
+    !which is available from https://atomsk.univ-lille.fr/. They have been adapted to fit cacmb.
+    subroutine disloop
+        !This subroutine applies the displacement field for a dislocation loop to all atoms and elements.
+        integer :: a1, a2, a3 !New directions
+        integer :: i, j, inod, ibasis, Npoints 
+        real(kind = dp) :: perimeter, angle, theta, omega, xA(3), xB(3), xC(3), u(3)
+        real(kind=dp), dimension(:,:), allocatable :: xloop !coordinate of points forming loop
+
+        if(allocated(xLoop)) deallocate(xLoop)
+
+        !Define new directions
+        select case(loop_normal)
+        case('x','X')
+            !Loop in (y,z) plane
+            a1 = 2
+            a2 = 3
+            a3 = 1
+        case('y','Y')
+            !Loop in (x,z) plane
+            a1 = 3
+            a2 = 1
+            a3 = 2
+        case default
+            !Loop in (x,y) plane
+            a1 = 1
+            a2 = 2
+            a3 = 3
+        end select
+
+        !Calculate loop perimeter
+        perimeter = 2.0_dp*pi*loop_radius
+
+        !Define the number of points forming the loop 
+        ! The following criteria are used as a trade-off between
+        ! good accuracy and computational efficiency:
+        ! - each dislocation segment should have a length of 5 angströms;
+        ! - the loop should contain at least 3 points
+        !    (for very small loops, this will result in segments shorter than 5 A);
+        ! - there should not be more than 100 points
+        !    (for very large loops, this will result in segments longer than 5 A).
+        Npoints = MAX( 3 , MIN( NINT(perimeter/5.d0) , 100 ) )
+
+        !angle between two consecutive points
+        theta = 2.0_dp*pi / dble(Npoints)
+
+        !allocate xLoop
+        allocate(xLoop(Npoints,3))
+        xLoop(:,:) = 0.0_dp
+
+        !Calculate the position of each point in the loop
+        angle = 0.0_dp
+        do i = 1, size(xLoop,1)
+            xLoop(i,a1) = centroid(a1) + loop_radius*dcos(angle)
+            xLoop(i,a2) = centroid(a2) + loop_radius*dsin(angle)
+            xLoop(i,a3) = centroid(a3)
+            ! Increment angle for next point
+            angle = angle + theta
+        end do
+
+        !Now actually calculate the displacement created by a loop for every atom
+        do i = 1, atom_num
+            u = 0.0_dp
+            omega = 0.0_dp
+            xC(:) = centroid - r_atom(:,i)
+
+            !Loop over all dislocation segments
+            do j = 1, size(xLoop,1)
+                !Coordinates of point A
+                if (j ==1) then
+                    xA(:) = xLoop(size(xLoop,1),:) - r_atom(:,i)
+                else
+                    xA(:) = xLoop(j-1, :) - r_atom(:,i)
+                end if
+
+                !Coordinates of point B
+                xB(:) = xLoop(j,:) - r_atom(:,i)
+
+                !Displacement due to solid angle
+                omega = omega + SolidAngle(xA, xB, xC)
+                !Displacement due to elasticity
+                u(:) = u(:) + DisloSeg_displacement_iso(xA, xB, burgers(:), poisson)
+            end do
+
+            !Total displacement
+            u=u(:) + burgers(:) * omega
+
+            !Apply displacement to atom
+            r_atom(:,i) = r_atom(:,i) + u(:)
+        end do
+
+        !Repeat for element nodes
+        do i = 1, ele_num
+            do inod =1, ng_node(lat_ele(i))
+                do ibasis = 1, basisnum(lat_ele(i))
+                    u = 0.0_dp
+                    omega = 0.0_dp
+                    xC(:) = centroid - r_node(:,ibasis,inod,i)
+
+                    !Loop over all dislocation segments
+                    do j = 1, size(xLoop,1)
+                        !Coordinates of point A
+                        if (j ==1) then
+                            xA(:) = xLoop(size(xLoop,1),:) - r_node(:,ibasis,inod,i)
+                        else
+                            xA(:) = xLoop(j-1, :) - r_node(:,ibasis,inod,i)
+                        end if
+
+                        !Coordinates of point B
+                        xB(:) = xLoop(j,:) - r_node(:,ibasis,inod,i)
+
+                        !Displacement due to solid angle
+                        omega = omega + SolidAngle(xA, xB, xC)
+                        !Displacement due to elasticity
+                        u(:) = u(:) + DisloSeg_displacement_iso(xA, xB, burgers(:), poisson)
+                    end do
+
+                    !Total displacement
+                    u=u(:) + burgers(:) * omega
+
+                    !Apply displacement to atom
+                    r_node(:,ibasis,inod,i) = r_node(:,ibasis,inod,i) + u(:)
+                end do
+            end do 
+        end do
+        return
+    end subroutine
+
+    !********************************************************
+    ! SOLIDANGLE
+    ! Calculate solid angle (normalized by 4*pi) used to obtain the 
+    ! displacement created by a dislocation triangle loop ABC at the origin
+    ! (field point = origin)
+    ! Ref.: Van Oosterom, A. and Strackee, J., 
+    !       The Solid Angle of a Plane Triangle, 
+    !       IEEE Transactions on Biomedical Engineering BME-30, 125 (1983).
+    !********************************************************
+    FUNCTION SolidAngle(xA, xB, xC) RESULT(Omega)
+    !
+    IMPLICIT NONE
+    !
+    ! Extremities of the triangle loop
+    REAL(dp),DIMENSION(3),INTENT(IN) :: xA, xB, xC
+    !
+    ! Solid angle (normalized by 4*pi)
+    REAL(dp):: omega
+    REAL(dp),PARAMETER:: factor=1.d0/(2.d0*pi)
+    !
+    REAL(dp) :: rA, rB, rC, numerator, denominator
+    !
+    rA = norm2(xA)
+    rB = norm2(xB)
+    rC = norm2(xC)
+    !
+    numerator = TRIPLE_PRODUCT( xA, xB, xC )
+    denominator = rA*rB*rC + DOT_PRODUCT( xA, xB )*rC &
+                & + DOT_PRODUCT( xB, xC )*rA + DOT_PRODUCT( xC, xA )*rB
+    !
+    omega = factor*ATAN2( numerator, denominator )
+    !
+    END FUNCTION SolidAngle
+
+    !********************************************************
+    ! DISLOSEG_DISPLACEMENT_ISO
+    ! Calculate displacement created by dislocation segment AB
+    ! once the solid angle part has been removed
+    ! Isotropic elastic calculation with nu Poisson coef.
+    ! Ref.: Eq. (1) in Barnett, D. M. 
+    !       The Displacement Field of a Triangular Dislocation Loop
+    !       Philos. Mag. A, 1985, 51, 383-387 
+    !********************************************************
+    FUNCTION DisloSeg_displacement_iso(xA, xB, b, nu) RESULT(u)
+    !
+    IMPLICIT NONE
+    REAL(dp),DIMENSION(3),INTENT(IN):: xA, xB ! Extremities of the segment
+    REAL(dp),DIMENSION(3),INTENT(IN):: b  ! Burgers vector
+    REAL(dp),INTENT(IN):: nu    ! Poisson coefficient
+    REAL(dp),DIMENSION(3):: u ! Displacement
+    REAL(dp):: rA, rB
+    REAL(dp),DIMENSION(1:3)::  tAB, nAB 
+    !
+    rA = norm2(xA)
+    rB = norm2(xB)
+    !
+    ! Tangent vector
+    tAB(:) = xB(:) - xA(:)
+    tAB(:) = tAB(:)/norm2(tAB)
+    !
+    ! Normal vector
+    nAB(:) = CROSS_PRODUCT(xA,xB)
+    nAB(:) = nAB(:)/norm2(nAB)
+    !
+    u(:) = ( -(1.d0-2.d0*nu)*CROSS_PRODUCT(b, tAB)* &
+         &  LOG( (rB + DOT_PRODUCT(xB,tAB))/(rA + DOT_PRODUCT(xA,tAB)) ) &
+         &  + DOT_PRODUCT(b,nAB)*CROSS_PRODUCT(xB/rB-xA/rA,nAB) ) &
+         &  /(8.d0*pi*(1.d0-nu))
+    !
+    END FUNCTION DisloSeg_displacement_iso
+!
+    !This code simply creates a planar vacancy cluster and does not apply the dislocation loop displacement field.
+    ! subroutine disloop
+    !     !This subroutine actually creates the dislocation loop.
+
+    !     real(kind=dp) :: neighbor_dis(loop_size), temp_box(6), dis
+    !     integer :: i, j, index(loop_size)
+
+    !     neighbor_dis(:) = HUGE(1.0_dp)
+    !     index(:) = 0
+
+    !     !First find the nearest atom to the centroid 
+    !     do i = 1, atom_num
+    !         if(norm2(r_atom(:,i) - centroid) < neighbor_dis(1)) then 
+    !             neighbor_dis(1) = norm2(r_atom(:,i) - centroid)
+    !             index(1) = i
+    !         end if
+    !     end do 
+
+    !     !Define a new box, this box tries to isolate all atoms on the plane of the atom 
+    !     !closest to the user defined centroid.
+    !     temp_box(:) = box_bd(:) 
+    !     temp_box(2*loop_normal) = r_atom(loop_normal,index(1)) + 10.0_dp**(-2.0_dp)
+    !     temp_box(2*loop_normal-1) = r_atom(loop_normal,index(1)) - 10.0_dp**(-2.0_dp)
+
+    !     !Now reset the list for the scanning algorithm
+    !     index(1) = 0
+    !     neighbor_dis(1) = HUGE(1.0_dp)
+
+    !     !Now scan over all atoms again and find the closest loop_size number of atoms to the initial atom 
+    !     !that reside on the same plane.
+
+    !     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)
+    !             do j = 1, loop_size
+    !                 !Check to see if it is closer than other atoms
+    !                 if (dis < neighbor_dis(j)) then 
+    !                     !Move values in the neighbor array and add the new neighbor
+    !                     if(j < loop_size) then 
+    !                         neighbor_dis(j+1:loop_size) = neighbor_dis(j:loop_size-1)
+    !                         index(j+1:loop_size) = index(j:loop_size-1)
+    !                     end if
+    !                     neighbor_dis(j) = dis
+    !                     index(j) = i
+    !                     exit
+    !                 end if
+    !             end do 
+    !         end if
+    !     end do 
+
+    !     !Now delete the atoms
+    !     call delete_atoms(loop_size, index)
+
+    !     return
+    ! end subroutine disloop
+
+end module opt_disl
\ No newline at end of file

From ba28734fd7ccf34b2a816910528c57d225cabcc4 Mon Sep 17 00:00:00 2001
From: Alex Selimov <alsel@tuta.io>
Date: Mon, 13 Jan 2020 16:54:52 -0500
Subject: [PATCH 14/17] Working changes to write lammps cac

---
 src/mode_create.f90 | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/mode_create.f90 b/src/mode_create.f90
index 6c465a3..b65d62e 100644
--- a/src/mode_create.f90
+++ b/src/mode_create.f90
@@ -306,9 +306,9 @@ module mode_create
              do inod = 1, 8
                 do i = 1,3
                     if(is_equal(cubic_cell(i, inod),0.0_dp)) then 
-                        adjustVar(i,inod) = -0.25_dp
+                        adjustVar(i,inod) = -0.5_dp
                     else 
-                        adjustVar(i, inod) = 0.25_dp
+                        adjustVar(i, inod) = 0.5_dp
                     end if
                 end do
             end do

From fbae3a15946d52fa55513ee2759ba18389913bc4 Mon Sep 17 00:00:00 2001
From: Alex Selimov <alsel@tuta.io>
Date: Mon, 13 Jan 2020 20:12:46 -0500
Subject: [PATCH 15/17] Working version of CACmb capable of building simple
 models for lammpsCAC

---
 src/Makefile        |  9 +++++----
 src/elements.f90    | 20 +++++++++++++++++-
 src/io.f90          | 49 +--------------------------------------------
 src/main.f90        | 11 ++++++++++
 src/mode_create.f90 | 18 +++--------------
 src/parameters.f90  |  2 ++
 6 files changed, 41 insertions(+), 68 deletions(-)

diff --git a/src/Makefile b/src/Makefile
index c84186e..e1a1145 100644
--- a/src/Makefile
+++ b/src/Makefile
@@ -1,6 +1,7 @@
+ 
 FC=ifort
-#FFLAGS=-mcmodel=large -g -O0 -stand f08 -fpe0 -traceback -check bounds,uninit -warn all -implicitnone -no-wrap-margin
-FFLAGS=-mcmodel=large -Ofast
+FFLAGS=-mcmodel=large -g -O0 -stand f08 -fpe0 -traceback -check bounds,uninit -warn all -implicitnone -no-wrap-margin
+#FFLAGS=-mcmodel=large -Ofast
 MODES=mode_create.o mode_merge.o mode_convert.o
 OBJECTS=main.o elements.o io.o subroutines.o functions.o atoms.o call_mode.o box.o $(MODES)
 
@@ -8,7 +9,7 @@ OBJECTS=main.o elements.o io.o subroutines.o functions.o atoms.o call_mode.o box
 .SUFFIXES: .c .f .f90 .F90 .o
 
 cacmb: $(OBJECTS) 
-	$(FC) $(FFLAGS) $(OBJECTS) -o $@
+	$(FC) $(FFLAGS) $(OBJECTS) parameters.o -o $@
 
 .f90.o: 
 	$(FC) $(FFLAGS) -c $< 
@@ -30,4 +31,4 @@ main.o io.o build_subroutines.o: elements.o
 call_mode.o : $(MODES)
 $(MODES) io.o: atoms.o box.o
 $(MODES) main.o : io.o
-testfuncs.o elements.o mode_create.o: subroutines.o
+testfuncs.o elements.o mode_create.o: subroutines.o
\ No newline at end of file
diff --git a/src/elements.f90 b/src/elements.f90
index 22e6eff..d6545d6 100644
--- a/src/elements.f90
+++ b/src/elements.f90
@@ -101,10 +101,28 @@ module elements
 
         real(kind=dp), dimension(3,max_ng_node), intent(out) :: cell_mat
 
+        integer :: inod, i
+        real(kind=dp), dimension(3,max_ng_node) :: adjustVar
+
+        adjustVar(:,:) = 0.0_dp
+
         select case(trim(ele_type))
 
         case('fcc')
-            cell_mat(:,1:8) = lapa * ((esize-1)*matmul(orient_mat, fcc_cell))
+            if(lmpcac) then 
+                do inod = 1, 8
+                    do i = 1,3
+                        if(is_equal(cubic_cell(i, inod),0.0_dp)) then 
+                            adjustVar(i,inod) = -0.5_dp
+                        else 
+                            adjustVar(i, inod) = 0.5_dp
+                        end if
+                    end do
+                end do
+                adjustVar(:,1:8) = matmul(fcc_mat, adjustVar(:,1:8))
+            end if
+            cell_mat(:, 1:8) = fcc_cell + adjustVar(:,1:8)
+            cell_mat(:,1:8) = lapa * ((esize-1)*matmul(orient_mat, cell_mat(:,1:8)))
         case default
             print *, "Element type ", trim(ele_type), " currently not accepted"
             stop
diff --git a/src/io.f90 b/src/io.f90
index dd49096..9d5537a 100644
--- a/src/io.f90
+++ b/src/io.f90
@@ -9,7 +9,7 @@ module io
 
     integer :: outfilenum = 0, infilenum = 0
     character(len=100) :: outfiles(10), infiles(10)
-    logical lmpcac
+    
 
    public 
    contains
@@ -78,53 +78,6 @@ module io
 
     end subroutine get_out_file 
 
-    subroutine get_in_file(filename)
-
-        implicit none
-
-        character(len=100), intent(in) :: filename
-        character(len=100) :: temp_infile
-        logical :: file_exists
-
-        !If no filename is provided then this function is called with none and prompts user input
-        if (filename=='none') then 
-            print *, "Please specify a filename with extension to read in:"
-            read(*,*) temp_infile
-        else
-            temp_infile = filename
-        end if
-
-        !Infinite loop which only exists if user provides valid filetype
-        do while(.true.)
-
-            !Check to see if file exists, if it doesn't then ask the user for another input
-            inquire(file=trim(temp_infile), exist=file_exists)
-            if (.not.file_exists) then 
-                print *, "The file ", temp_infile, "does not exist, please input an existing file to read in."
-                read(*,*) temp_infile
-                cycle
-            end if
-
-            if (scan(temp_infile,'.',.true.) == 0) then 
-                print *, "No extension included on filename, please type a full filename that includes an extension."
-                read(*,*) temp_infile
-                cycle
-            end if
-            select case(temp_infile(scan(temp_infile,'.',.true.)+1:))
-            case('cac')
-                infilenum=infilenum+1
-                infiles(infilenum) =  temp_infile
-                exit
-            case default    
-                print *, "File type: ", trim(temp_infile(scan(temp_infile,'.',.true.):)), " not currently accepted. ", &
-                         "please input a filename with extension from following list: cac."
-                read(*,*) temp_infile
-
-            end select
-        end do
-
-    end subroutine get_in_file 
-
     subroutine write_out
         !This subroutine loops over alll of the outfile types defined and calls the correct writing subroutine
 
diff --git a/src/main.f90 b/src/main.f90
index a681d33..f90b731 100644
--- a/src/main.f90
+++ b/src/main.f90
@@ -27,6 +27,17 @@ program main
     ! Command line parsing
     arg_num = command_argument_count()
 
+    !First check if we are writing out to lammpscac format by looping over all arguments
+    do i = 1, arg_num
+        call get_command_argument(i, argument)
+        select case(argument(scan(argument,'.',.true.)+1:))
+        case('cac')
+            lmpcac = .true.
+        case default
+            continue
+        end select
+    end do 
+
     !Determine if a mode is being used and what it is. The first argument has to be the mode
     !if a mode is being used
     call get_command_argument(1, argument)
diff --git a/src/mode_create.f90 b/src/mode_create.f90
index d5995e8..a669f93 100644
--- a/src/mode_create.f90
+++ b/src/mode_create.f90
@@ -91,20 +91,6 @@ module mode_create
             end do 
         else
 
-            if(lmpcac) then 
-                 do inod = 1, 8
-                    do i = 1,3
-                        if(is_equal(cubic_cell(i, inod),0.0_dp)) then 
-                            adjustVar(i,inod) = -0.5_dp
-                        else 
-                            adjustVar(i, inod) = 0.5_dp
-                        end if
-                    end do
-                end do
-            else
-                adjustVar(:,:)=0.0_dp
-            end if
-
             call cell_init(lattice_parameter, esize, element_type, orient, cell_mat) 
             !If the user doesn't pass any build instructions than we just put the cell mat into the element_array 
             call alloc_ele_arrays(1,0)
@@ -333,6 +319,8 @@ module mode_create
                     end if
                 end do
             end do
+
+            adjustVar(:,1:8) = matmul(orient,matmul(fcc_mat,adjustVar(:,1:8)))
         else
             adjustVar(:,:)=0.0_dp
         end if
@@ -465,7 +453,7 @@ module mode_create
                         end do
 
                         if(all(node_in_bd)) then 
-                            lat_ele_num = lat_ele_num+1
+                            lat_ele_num = lat_ele_num+1 
                             do inod = 1, 8
                                 r_lat(:,inod,lat_ele_num) = temp_nodes(:,1,inod) + adjustVar(:,inod)
                             end do
diff --git a/src/parameters.f90 b/src/parameters.f90
index 0443622..9e9d9d8 100644
--- a/src/parameters.f90
+++ b/src/parameters.f90
@@ -5,4 +5,6 @@ module parameters
     integer, parameter :: dp= selected_real_kind(15,307)
     real(kind=dp), parameter :: lim_zero = epsilon(1.0_dp), &
                                 lim_large = huge(1.0_dp)
+    logical, save :: lmpcac
+    
 end module parameters
\ No newline at end of file

From 55fbe679e5ef89ef516014d9d97b8d71dd4a5ff4 Mon Sep 17 00:00:00 2001
From: Alex Selimov <alsel@tuta.io>
Date: Wed, 15 Jan 2020 08:12:48 -0500
Subject: [PATCH 16/17] Quick fix to argument positioning

---
 src/mode_merge.f90 | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/mode_merge.f90 b/src/mode_merge.f90
index ca9d068..cfb4273 100644
--- a/src/mode_merge.f90
+++ b/src/mode_merge.f90
@@ -77,15 +77,15 @@ module mode_merge
         end do
 
         !Set argpos accordingly
-        arg_pos = 3+in_num+1
+        arg_pos = 3+in_num
         !Now options loop
         !Check for optional keywords
         do while(.true.)
             if(arg_pos > command_argument_count()) exit
             !Pull out the next argument which should either be a keyword or an option
+            arg_pos=arg_pos+1
             call get_command_argument(arg_pos, textholder)
             textholder=adjustl(textholder)
-            arg_pos=arg_pos+1
 
             !Choose what to based on what the option string is 
             select case(trim(textholder))

From 849da1d24ac1767bd72757261364039eea26c5aa Mon Sep 17 00:00:00 2001
From: Alex Selimov <alsel@tuta.io>
Date: Wed, 15 Jan 2020 09:39:30 -0500
Subject: [PATCH 17/17] Changes to how the adjustment to nodal positions is
 performed for lammpscac output

---
 src/elements.f90    | 19 +++++++++++-----
 src/io.f90          | 38 ++++++++++++++++++++++++-------
 src/main.f90        | 11 ---------
 src/mode_create.f90 | 55 +++++++++++++++++++++------------------------
 4 files changed, 69 insertions(+), 54 deletions(-)

diff --git a/src/elements.f90 b/src/elements.f90
index 6b41e57..b4e59be 100644
--- a/src/elements.f90
+++ b/src/elements.f90
@@ -9,7 +9,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(:) !Element siz
+    integer, allocatable :: size_ele(:), lat_ele(:), sbox_ele(:) !Element size
     real(kind=dp), allocatable :: r_node(:,:,:,:) !Nodal position array
 
     integer, save :: ele_num !Number of elements
@@ -39,6 +39,7 @@ module elements
     !This can be easily increased with no change to efficiency
     integer :: max_basisnum, basisnum(10) !Max basis atom number, number of basis atoms in each lattice type 
     integer :: basis_type(10,10)
+    real(kind=dp) :: lapa(10)
 
     public 
     contains
@@ -140,7 +141,7 @@ module elements
         
         !Allocate element arrays
         if(n > 0) then 
-            allocate(type_ele(n), size_ele(n), lat_ele(n), r_node(3,max_basisnum, max_ng_node,n), &
+            allocate(type_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
@@ -179,12 +180,17 @@ module elements
             allocate(temp_int(n+ele_num+buffer_size))
             temp_int(1:ele_size) = lat_ele
             temp_int(ele_size+1:) = 0
-            call move_alloc(temp_int(1:ele_size), lat_ele)
+            call move_alloc(temp_int, lat_ele)
 
             allocate(temp_int(n+ele_num+buffer_size))
             temp_int(1:ele_size) = size_ele
             temp_int(ele_size+1:) = 0
-            call move_alloc(temp_int(1:ele_size), size_ele)
+            call move_alloc(temp_int, size_ele)
+
+            allocate(temp_int(n+ele_num+buffer_size))
+            temp_int(1:ele_size) = lat_ele
+            temp_int(ele_size+1:) = 0
+            call move_alloc(temp_int, sbox_ele)
 
             allocate(char_temp(n+ele_num+buffer_size))
             char_temp(1:ele_size) = type_ele
@@ -210,9 +216,9 @@ module elements
         end if
     end subroutine
 
-    subroutine add_element(type, size, lat, r)
+    subroutine add_element(type, size, lat, sbox, r)
         !Subroutine which adds an element to the element arrays
-        integer, intent(in) :: size, lat
+        integer, intent(in) :: size, lat, sbox
         character(len=100), intent(in) :: type
         real(kind=dp), intent(in) :: r(3, max_basisnum, max_ng_node)
 
@@ -220,6 +226,7 @@ module elements
         type_ele(ele_num) = type
         size_ele(ele_num) = size
         lat_ele(ele_num) = lat
+        sbox_ele(ele_num) = sbox
        r_node(:,:,:,ele_num) = r(:,:,:)
        node_num = node_num + ng_node(lat)
 
diff --git a/src/io.f90 b/src/io.f90
index 4f9b74c..e5b098d 100644
--- a/src/io.f90
+++ b/src/io.f90
@@ -204,7 +204,7 @@ module io
         !This subroutine writes out a .lmp style dump file
         character(len=100), intent(in) :: file
         integer :: write_num, i, inod, ibasis
-        real(kind=dp) :: mass
+        real(kind=dp) :: mass, fcc_adjust(3,8), local_adjust(3,8), rout(3)
 
 1       format(i16, '   Eight_Node', 4i16)
 2       format(i16, '   Atom', 4i16)
@@ -242,16 +242,32 @@ module io
         write(11, '(a)') 'CAC Elements'
         write(11, '(a)') ' '
 
+        !Set up the nodal adjustment variables for all the different element types. This adjusts the node centers
+        !from the center of the unit cell (as formulated in this code) to the corners of the unit cells
+         do inod = 1, 8
+            do i = 1,3
+                if(is_equal(cubic_cell(i, inod),0.0_dp)) then 
+                    fcc_adjust(i,inod) = -0.5_dp
+                else 
+                    fcc_adjust(i, inod) = 0.5_dp
+                end if
+            end do
+        end do
+        fcc_adjust = matmul(fcc_mat, fcc_adjust)
+
         !Write element nodal positions
         do i = 1, ele_num
             select case(trim(adjustl(type_ele(i))))
             case('fcc')
+                !Now orient the current adjustment vector to the correct orientation
+                local_adjust = matmul(sub_box_ori(:,:,sbox_ele(i)), fcc_adjust) * lapa(lat_ele(i))
                 !The first entry is the element specifier 
                 write(11,1) i, basisnum(lat_ele(i)), size_ele(i), size_ele(i), size_ele(i)
                 do ibasis = 1, basisnum(lat_ele(i))
                     do inod = 1, 8
                         !Nodal information for every node
-                        write(11,3) inod, ibasis, basis_type(ibasis,lat_ele(i)), r_node(:,ibasis,inod,i)
+                        rout = r_node(:,ibasis,inod,i) + local_adjust(:,inod)
+                        write(11,3) inod, ibasis, basis_type(ibasis,lat_ele(i)), rout
                     end do
                 end do
             end select
@@ -486,6 +502,8 @@ module io
         write(11,*) lattice_types, (basisnum(i), i = 1, lattice_types), (ng_node(i), i = 1, lattice_types)
         !Now for every lattice type write the basis atom types
         write(11,*) ((basis_type(i,j), i = 1, basisnum(j)), j = 1, lattice_types)
+        !Now for every lattice type write the lattice parameters
+        write(11,*) (lapa(i), i = 1, lattice_types)
 
         !Now write the numbers of elements and atoms
         write(11,*) atom_num, ele_num
@@ -498,7 +516,7 @@ module io
         !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), type_ele(i)
+            write(11, *) 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)
@@ -582,7 +600,7 @@ module io
         real(kind = dp), dimension(6), intent(out) :: temp_box_bd
 
         integer :: i, j, k, n, inod, ibasis, type, size, in_atoms, in_eles, new_atom_types, &
-                   new_type_to_type(10), new_lattice_types
+                   new_type_to_type(10), new_lattice_types, sbox
         character(len=100) :: etype
         real(kind=dp) ::  r(3), newdisplace(3)
         real(kind=dp), allocatable :: r_innode(:,:,:)
@@ -623,8 +641,6 @@ module io
         sub_box_array_bd(:,1,sub_box_num+1:) = sub_box_array_bd(:,1,sub_box_num+1:) + atom_num
         sub_box_array_bd(:,2,sub_box_num+1:) = sub_box_array_bd(:,2,sub_box_num+1:) + ele_num
 
-        sub_box_num = sub_box_num + n
-
         !Read in the number of atom types and all their names
         read(11, *) new_atom_types, (new_type_to_name(i), i = 1, new_atom_types)
         !Now fit these into the global list of atom types, after this new_type_to_type is the actual global 
@@ -646,6 +662,8 @@ module io
                 basis_type(i,j) = new_type_to_type(basis_type(i,j))
             end do 
         end do 
+        !Read the lattice parameters for every lattice type
+        read(11,*) (lapa(i), i = lattice_types+1, lattice_types+new_lattice_types)
         !Read number of elements and atoms and allocate arrays
         read(11, *) in_atoms, in_eles
         call grow_ele_arrays(in_eles, in_atoms)
@@ -659,7 +677,7 @@ module io
 
         !Read the elements
         do i = 1, in_eles
-            read(11, *) n, type, size, etype
+            read(11, *) n, type, size, sbox, etype
             do inod = 1, ng_node(type)
                 do ibasis =1, basisnum(type)
                     read(11,*) j, k, r_innode(:, ibasis, inod)
@@ -667,7 +685,7 @@ module io
                 end do
             end do
             type = type + lattice_types
-            call add_element(etype, size, type, r_innode)
+            call add_element(etype, size, type, sbox+n, r_innode)
         end do
 
         !Close the file being read
@@ -676,5 +694,9 @@ module io
         !Only increment the lattice types if there are elements, if there are no elements then we 
         !just overwrite the arrays
         if(in_eles > 0) lattice_types = lattice_types + new_lattice_types
+
+
+        sub_box_num = sub_box_num + n
+
     end subroutine read_mb
 end module io
diff --git a/src/main.f90 b/src/main.f90
index ccb4d31..ed34c53 100644
--- a/src/main.f90
+++ b/src/main.f90
@@ -28,17 +28,6 @@ program main
     ! Command line parsing
     arg_num = command_argument_count()
 
-    !First check if we are writing out to lammpscac format by looping over all arguments
-    do i = 1, arg_num
-        call get_command_argument(i, argument)
-        select case(argument(scan(argument,'.',.true.)+1:))
-        case('cac')
-            lmpcac = .true.
-        case default
-            continue
-        end select
-    end do 
-
     !Determine if a mode is being used and what it is. The first argument has to be the mode
     !if a mode is being used
     call get_command_argument(1, argument)
diff --git a/src/mode_create.f90 b/src/mode_create.f90
index 9b2cd73..9463e92 100644
--- a/src/mode_create.f90
+++ b/src/mode_create.f90
@@ -12,7 +12,7 @@ module mode_create
 
     character(len=100) :: name, element_type
     real(kind = dp) :: lattice_parameter, orient(3,3), cell_mat(3,8), box_len(3), basis(3,3), origin(3), maxlen(3), &
-                       orient_inv(3,3), box_vert(3,8), maxbd(3), lattice_space(3), duplicate(3), adjustVar(3,8)
+                       orient_inv(3,3), box_vert(3,8), maxbd(3), lattice_space(3), duplicate(3)
     integer :: esize, ix, iy, iz, box_lat_vert(3,8), lat_ele_num, lat_atom_num, bd_in_lat(6), &
                basis_pos(3,10)
     logical :: dup_flag, dim_flag
@@ -105,7 +105,7 @@ module mode_create
                 box_bd(2*i) = maxval(r_node_temp(i,:,:))
                 box_bd(2*i-1) = origin(i)
             end do
-            call add_element(element_type, esize, 1, r_node_temp)
+            call add_element(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
@@ -141,7 +141,7 @@ module mode_create
                             r_node_temp(:,ibasis,inod) = (r_lat(:,inod,i)*lattice_parameter)+basis_pos(:,ibasis)
                         end do
                     end do 
-                    call add_element(element_type, esize, 1, r_node_temp)
+                    call add_element(element_type, esize, 1, 1, r_node_temp)
                 end do 
             end if
         end if
@@ -199,7 +199,20 @@ module mode_create
                     call get_command_argument(arg_pos, orient_string, arglen)
                     if (arglen==0) STOP "Missing orientation in orient command of mode create"
                     arg_pos = arg_pos+1
-                    call parse_ori_vec(orient_string, orient(i,:))
+                    ori_pos=2
+                    do j = 1,3
+                        if (orient_string(ori_pos:ori_pos) == '-') then 
+                            ori_pos = ori_pos + 1
+                            read(orient_string(ori_pos:ori_pos), *, iostat=stat) orient(i,j)
+                            if (stat>0) STOP "Error reading orient value"
+                            orient(i,j) = -orient(i,j)
+                            ori_pos = ori_pos + 1
+                        else
+                            read(orient_string(ori_pos:ori_pos), *, iostat=stat) orient(i,j)
+                            if(stat>0) STOP "Error reading orient value"
+                            ori_pos=ori_pos + 1
+                        end if
+                    end do
                 end do
 
 
@@ -234,6 +247,7 @@ module mode_create
                 exit
             end select
         end do
+
         !Calculate the lattice periodicity length in lattice units
         do i = 1, 3
             lattice_space(i) = norm2(orient(i,:))
@@ -264,8 +278,9 @@ module mode_create
         !Now normalize the orientation matrix
         orient =  matrix_normal(orient,3)
 
-        !Set lattice_num to 1 
+        !Set lattice_num to 1 and add the lattice_parameter to the elements module lattice paramter variable
         lattice_types = 1
+        lapa(1) = lattice_parameter
 
         !If we haven't defined a basis then define the basis and add the default basis atom type and position
         if (basisnum(1) == 0) then
@@ -286,7 +301,7 @@ module mode_create
         !Internal variables
         integer :: i,  inod, bd_in_lat(6), bd_in_array(6), ix, iy, iz,  numlatpoints, ele(3,8), rzero(3), &
                    vlat(3), temp_lat(3,8), m, n, o
-        real(kind=dp) :: v(3), temp_nodes(3,1,8), adjustVar(3,8)
+        real(kind=dp) :: v(3), temp_nodes(3,1,8)
         real(kind=dp), allocatable :: resize_lat_array(:,:)
         logical, allocatable :: lat_points(:,:,:)
         logical :: node_in_bd(8)
@@ -294,23 +309,6 @@ module mode_create
         !Do some value initialization
         max_esize = esize
 
-        !If we are writing out to lammpscac format then we have to adjust the nodal positions
-
-        if(lmpcac) then 
-             do inod = 1, 8
-                do i = 1,3
-                    if(is_equal(cubic_cell(i, inod),0.0_dp)) then 
-                        adjustVar(i,inod) = -0.5_dp
-                    else 
-                        adjustVar(i, inod) = 0.5_dp
-                    end if
-                end do
-            end do
-
-            adjustVar(:,1:8) = matmul(orient,matmul(fcc_mat,adjustVar(:,1:8)))
-        else
-            adjustVar(:,:)=0.0_dp
-        end if
         !First find the bounding lattice points (min and max points for the box in each dimension)
         numlatpoints = 1
         do i = 1, 3
@@ -328,6 +326,7 @@ module mode_create
             continue
         end select
 
+
         !Loop over all of lattice points within the boundary, we choose between two loops. One for the atomistic case 
         !and one for the regular case
         if (esize==2) then  
@@ -440,11 +439,9 @@ module mode_create
                         end do
 
                         if(all(node_in_bd)) then 
-                            lat_ele_num = lat_ele_num+1 
-                            do inod = 1, 8
-                                r_lat(:,inod,lat_ele_num) = temp_nodes(:,1,inod) + adjustVar(:,inod)
-                            end do
-                            
+                            lat_ele_num = lat_ele_num+1
+                            r_lat(:,:,lat_ele_num) = temp_nodes(:,1,:)
+
                             !Now set all the lattice points contained within an element to false
                             do o = minval(temp_lat(3,:)), maxval(temp_lat(3,:))
                                 do n = minval(temp_lat(2,:)), maxval(temp_lat(2,:))
@@ -510,4 +507,4 @@ module mode_create
     end subroutine error_message
 
 
-end module mode_create
+end module mode_create
\ No newline at end of file