First working version of reading in .cac format. Works for two separate orientations

src/call_option.f90

@@ -34,7 +34,7 @@ subroutine call_option(option, arg_pos)
     case('-delete')
         call run_delete(arg_pos)
     case('-set_cac')
-        arg_pos = arg_pos+3
+        arg_pos=arg_pos +3
     case default
         print *, 'Option ', trim(adjustl(option)), ' is not currently accepted.'
         stop 3 

src/elements.f90

@@ -646,7 +646,7 @@ module elements
             if ((basisnum(j) == in_bnum).and.(ng_node(j) == in_ngnodes).and.(is_equal(lapa(j),in_lapa))) then 
                 !Now check lattice level variables 
                 do ibasis = 1, basisnum(j)
-                    if(basis_type(ibasis,j) /= in_btypes(ibasis)) cycle old_loop
+                    if(basis_type(ibasis,j) /= in_btypes(ibasis)) cycle lat_loop
                 end do
                 lat_type = j
                 exit lat_loop
@@ -657,9 +657,10 @@ module elements
         if( lat_type == 0) then 
             lattice_types = lattice_types + 1
             basisnum(lattice_types) = in_bnum
-            basis_types(:,lattice_types) = in_btypes
+            basis_type(:,lattice_types) = in_btypes
             ng_node(lattice_types) = in_ngnodes
             lapa(lattice_types) = in_lapa
+            lat_type = lattice_types
         end if
 
     end subroutine lattice_map

src/io.f90

@@ -591,7 +591,7 @@ module io
             end if
 
             select case(temp_infile(scan(temp_infile,'.',.true.)+1:))
-            case('restart', 'mb')
+            case('restart', 'mb', 'cac')
                 infilenum=infilenum+1
                 infiles(infilenum) =  temp_infile
                 exit
@@ -618,6 +618,8 @@ module io
             call read_mb(infiles(i), displace, temp_box_bd)
         case('restart')
             call read_pycac(infiles(i), displace, temp_box_bd)
+        case('cac')
+            call read_lmpcac(infiles(i), displace, temp_box_bd)
         case default
             print *, "The extension ", trim(adjustl(outfiles(i)(scan(outfiles(i),'.',.true.)+1:))), &
                      " is not accepted for writing. Please select from: mb and try again"
@@ -956,10 +958,11 @@ module io
         real(kind=dp), dimension(3), intent(in) :: displace
         real(kind = dp), dimension(6), intent(out) :: temp_box_bd
         
-        character(len=100) :: textholder, element_type, esize
+        character(len=100) :: textholder, element_type
         character(len=2) :: atom_species
-        integer :: i, j, k, ele_in, type_in, type_map(10), in_basis, node_types(10,8), inod, ibasis, in_basis_types(10)
+        integer :: i, j, k, ele_in, type_in, type_map(10), in_basis, node_types(10,8), inod, ibasis, in_basis_types(10), esize, &
-        real(kind=dp) :: mass, r_in(3,10,8), lat_vec(3,3), in_ori(3,3), temp_box_bd(6), newdisplace(3)
+                   lat_type  
+        real(kind=dp) :: mass, r_in(3,10,8), lat_vec(3,3), in_ori(3,3), newdisplace(3)
         
         !First check to make sure that we have set the needed variables
         if(is_equal(in_lapa,0.0_dp).or.(in_lattice_type=='')) then 
@@ -969,22 +972,21 @@ module io
         !Open the file
         open(unit=11, file=trim(adjustl(file)), action='read',position='rewind')
         
-        !Now initialiaze some important variables
+        !Now initialiaze some important variables if they aren't defined
-        max_basis_num = 10
+        if (max_basisnum==0) max_basisnum = 10
+        if (max_ng_node==0) max_ng_node=8
 
         !Read header information
         read(11, *) textholder
-        read(11, *) textholder
 
         !Read number of elements
         read(11, *) ele_in, textholder 
         read(11, *) type_in, textholder
 
         !Read box_boundaries
-        read(11,*) textholder
+        read(11,*) temp_box_bd(1:2), textholder
-        read(11,*) temp_box_bd(1:2), texholder
+        read(11,*) temp_box_bd(3:4), textholder
-        read(11,*) temp_box_bd(3:4), texholder
+        read(11,*) temp_box_bd(5:6), textholder 
-        read(11,*) temp_box_bd(5:6), texholder 
 
         !Shift the box boundaries if needed
         do i = 1, 3
@@ -1016,7 +1018,6 @@ module io
 
         !Read useless information
         read(11,*) textholder 
-        read(11,*) textholder 
 
         !Read atomic masses
         do i = 1, type_in
@@ -1026,13 +1027,11 @@ module io
         end do
 
         !Read useless info
-        do i = 1, 3
         read(11,*) textholder
-        end do
 
         !Start the reading loop
         do i = 1, ele_in
-            read(11,*) j, ele, element_type, in_basis, esize
+            read(11,*) j, element_type, in_basis, esize
             !Check to see if we need to grow the max_basis_num
             select case(trim(adjustl(element_type)))
             case('Eight_Node')
@@ -1044,9 +1043,9 @@ module io
                 !Now calculate the lattice vectors and shift the nodal points from the corners to the center of the unit cell
                 !Please check the nodal numbering figure in the readme in order to understand which nodes are used for the
                 !calculation
-                lat_vec(:,1) = (r_node(:,1,2) - r_node(:,1,1))/esize
+                lat_vec(:,1) = (r_in(:,1,2) - r_in(:,1,1))/(2*esize)
-                lat_vec(:,2) = (r_node(:,1,4) - r_node(:,1,1))/esize
+                lat_vec(:,2) = (r_in(:,1,4) - r_in(:,1,1))/(2*esize)
-                lat_vec(:,3) = (r_node(:,1,5) - r_node(:,1,1))/esize
+                lat_vec(:,3) = (r_in(:,1,5) - r_in(:,1,1))/(2*esize)
                 
                 !Now shift all the nodal positions 
                 select case(trim(adjustl(in_lattice_type)))
@@ -1068,7 +1067,7 @@ module io
                 call lattice_map(in_basis, in_basis_types, 8, in_lapa, lat_type)
 
                 !Now add the element 
-                call add_element(in_lattice_type, esize, lat_type, sub_box_num, r_in)
+                call add_element(in_lattice_type, esize, lat_type, sub_box_num, r_in(:,1:max_basisnum,1:max_ng_node))
                 
             case('Atom')
                 read(11, *) inod, ibasis, in_basis_types(ibasis), r_in(:,1,1)
@@ -1083,16 +1082,23 @@ module io
         integer, intent(in) :: apos
         integer :: arglen, arg_pos
 
+        character(len=100) :: textholder
+        
         arg_pos = apos + 1 
-        call get_command_argument(arg_pos, in_lapa, arglen)
+        call get_command_argument(arg_pos, textholder, arglen)
         if (arglen==0) then 
             print *, "Missing lattice parameter for set_input_lat"
         end if
+        read(textholder,*) in_lapa
+        print *, in_lapa
 
         arg_pos = arg_pos + 1
-        call get_command_argument(arg_pos, in_lattice_type, arglen)
+        call get_command_argument(arg_pos, textholder, arglen)
         if (arglen==0) then 
             print *, "Missing lattice type for set_input_lat"
         end if
-    end subroutine set_input_lat(arg_pos)
+        read(textholder,*) in_lattice_type
+        print *, in_lattice_type
+
+    end subroutine set_cac
 end module io

src/main.f90

@@ -61,10 +61,8 @@ program main
         case('-wrap')
             wrap_flag=.true.
         
-        !This gives necessary information in order to correctly read .cac files
         case('-set_cac')
             call set_cac(i)
-
         end select 
     end do 
     !Determine if a mode is being used and what it is. The first argument has to be the mode