Latest working changes to read-CAC

src/elements.f90

@@ -294,8 +294,6 @@ module elements
 
         integer :: i
 
-        max_ng_node = 0
-
         do i=1, n
             select case(trim(adjustl(element_types(i))))
             case('fcc')

src/io.f90

@@ -776,7 +776,7 @@ module io
         integer :: i, inod, ibasis, j, k, in_eles, in_atoms, ele_types, in_lat_num, in_atom_types, &
                    atom_type_map(10), etype_map(10), etype, lat_type, new_lattice_map(10), &
                    atom_type
-        real(kind=dp) :: newdisplace(3), r_in(3,1,8), r_in_atom(3), new_displace(3)
+        real(kind=dp) :: newdisplace(3), r_in(3,1,8), r_in_atom(3)
         character(len=100) :: textholder, in_lattype_map(10) 
         character(len=2) :: atomic_element
         !First open the file
@@ -950,7 +950,7 @@ module io
         end if
     end subroutine read_pycac
 
-    subroutine read_cac(file, displace, temp_box_bd)
+    subroutine read_lmpcac(file, displace, temp_box_bd)
         !This subroutine is used to read .cac files which are used with the lammpsCAC format
         character(len=100), intent(in) :: file
         real(kind=dp), dimension(3), intent(in) :: displace
@@ -958,8 +958,8 @@ module io
         
         character(len=100) :: textholder, element_type, esize
         character(len=2) :: atom_species
-        integer :: i, j, ele_in, type_in, type_map(10), in_basis, node_types(10,8), inod, ibasis, in_basis_types(10)
+        integer :: i, j, k, ele_in, type_in, type_map(10), in_basis, node_types(10,8), inod, ibasis, in_basis_types(10)
-        real(kind=dp) :: mass, r_in(3,10,8), lat_vec(3), in_ori(3,3)
+        real(kind=dp) :: mass, r_in(3,10,8), lat_vec(3,3), in_ori(3,3), temp_box_bd(6), 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,6 +969,9 @@ module io
         !Open the file
         open(unit=11, file=trim(adjustl(file)), action='read',position='rewind')
         
+        !Now initialiaze some important variables
+        max_basis_num = 10
+
         !Read header information
         read(11, *) textholder
         read(11, *) textholder
@@ -979,9 +982,38 @@ module io
 
         !Read box_boundaries
         read(11,*) textholder
-        read(11,*) box_bd(1:2), texholder
+        read(11,*) temp_box_bd(1:2), texholder
-        read(11,*) box_bd(3:4), texholder
+        read(11,*) temp_box_bd(3:4), texholder
-        read(11,*) box_bd(5:6), texholder 
+        read(11,*) temp_box_bd(5:6), texholder 
+
+        !Shift the box boundaries if needed
+        do i = 1, 3
+            if (abs(displace(i)) > lim_zero) then 
+                newdisplace(i) = displace(i) - temp_box_bd(2*i-1)
+            else
+                newdisplace(i)=displace(i)
+            end if 
+            temp_box_bd(2*i-1) = temp_box_bd(2*i-1) + newdisplace(i)
+            temp_box_bd(2*i) = temp_box_bd(2*i) + newdisplace(i)
+        end do 
+
+        !Grow box boundaries
+        call grow_box(temp_box_bd)
+
+        !Allocate sub_box 
+        if (sub_box_num == 0) then 
+            call alloc_sub_box(1)
+        else
+            call grow_sub_box(1)
+        end if
+
+        !Because orientations and other needed sub_box information isn't really 
+        !present within the .cac file we just default a lot of it.
+        sub_box_ori(:,:,sub_box_num+1) = identity_mat(3)
+        sub_box_bd(:, sub_box_num+1) = temp_box_bd
+        sub_box_num = sub_box_num + 1
+
+
         !Read useless information
         read(11,*) textholder 
         read(11,*) textholder 
@@ -1001,25 +1033,50 @@ module io
         !Start the reading loop
         do i = 1, ele_in
             read(11,*) j, ele, 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')
-                    !Read in all the data
+                !Read in all the data
-                    do j = 1, 8*in_basis 
+                do j = 1, 8*in_basis 
-                        read(11, *) inod, ibasis, in_basis_types(ibasis), r_in(:,ibasis,inod)
+                    read(11, *) inod, ibasis, in_basis_types(ibasis), r_in(:,ibasis,inod)
-                    end do
+                end do
-                     
-                    !Now calculate the orientation matrix based on the lattice type.
-                    lat_vec = r_in(:,1,2) - r_in(:,1,1) 
-                    lat_vec = lat_vec / norm2(lat_vec)
-                    
-                    !Now figure out if is an existing lattice_type
-                    call lattice_map(in_basis, in_basis_types, 8, in_lapa, lat_type)
-
                 
+                !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(:,2) = (r_node(:,1,4) - r_node(:,1,1))/esize
+                lat_vec(:,3) = (r_node(:,1,5) - r_node(:,1,1))/esize
+               
+                !Now shift all the nodal positions 
+                select case(trim(adjustl(in_lattice_type)))
+                case('fcc','FCC')
+                    do ibasis = 1, in_basis 
+                        r_in(:,ibasis,1) = r_in(:,ibasis,1) + lat_vec(:,1) + lat_vec(:,2) + lat_vec(:,3) + newdisplace
+                        r_in(:,ibasis,2) = r_in(:,ibasis,2) - lat_vec(:,1) + lat_vec(:,2) + lat_vec(:,3) + newdisplace
+                        r_in(:,ibasis,3) = r_in(:,ibasis,3) - lat_vec(:,1) - lat_vec(:,2) + lat_vec(:,3) + newdisplace
+                        r_in(:,ibasis,4) = r_in(:,ibasis,4) + lat_vec(:,1) - lat_vec(:,2) + lat_vec(:,3) + newdisplace
+                        r_in(:,ibasis,5) = r_in(:,ibasis,5) + lat_vec(:,1) + lat_vec(:,2) - lat_vec(:,3) + newdisplace
+                        r_in(:,ibasis,6) = r_in(:,ibasis,6) - lat_vec(:,1) + lat_vec(:,2) - lat_vec(:,3) + newdisplace
+                        r_in(:,ibasis,7) = r_in(:,ibasis,7) - lat_vec(:,1) - lat_vec(:,2) - lat_vec(:,3) + newdisplace
+                        r_in(:,ibasis,8) = r_in(:,ibasis,8) + lat_vec(:,1) - lat_vec(:,2) - lat_vec(:,3) + newdisplace
+                    end do
+                case default
+                    print *, in_lattice_type, " is not an accepted lattice type. Please select from: fcc"
+                end select 
+                !Now map it to either an existing or new lattice type
+                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)
+                
+            case('Atom')
+                read(11, *) inod, ibasis, in_basis_types(ibasis), r_in(:,1,1)
+                call add_atom(in_basis_types(ibasis), sub_box_num, r_in(:,1,1)) 
             end select
         end do
 
-    end subroutine read_cac
+    end subroutine read_lmpcac
 
     subroutine set_cac(apos)
         !This code parses input values